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FOREWORD 

This  report  incorporates  the  effects  of  shear  wave  from  the  elastic  ocean 
floor  to  the  normal  mode  transmission  loss  model.  The  possible  range 
dependence  of  the  acoust-ie  pror^o^ties  of  the  ocean  is  modeled  with  a  modified 
version  of  the  adiabatic  approximation  for  slow  variations.  This  newly 
developed  model  proves  itself  very  useful  for  performance  estimation  of 
acoustic  mines  and  sonar  systems  that  operate  at  low  frequencies.  This  work 
was  funded  by  the  Independent  Research  Board  of  the  Naval  Surface  Warfare 
Center . 
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CHAPTER  1 
INTRODUCTION 


An  electromagnetic  wave  in  vacuum  travels  a  distance  of  nearly  300,000 

kilometers  in  just  one  second.  At  the  same  time  sound  has  only  traveled  about 

340  meters.  However,  the  electromagnetic  wave  attenuates  very  rapidly  in 

water!  Also,  the  attenuation  coefficient  of  both  waves  varies  significantly 

with  frequency,  but  low-frequency  sound  is  the  best  form  of  underwater 

2 

radiation  known  to  man. 

An  acoustic  wave  incident  from  the  water  to  the  solid  sediments  of  the 
bottom  create  a  compressional  (longitudinal)  wave  and  a  shear  (transverse) 

3*6 

wave  as  shown  by  the  use  of  the  plane  wave  approximation  of  sound.  This,  so 
called,  birefraction  phenomena  is  also  observed  with  electromagnetic  waves  in 
calcite  and  it  is  described  by  a  medium  with  two  indices  of  refraction.  The 
same  phenomenon  applies  to  seismo-acoustics. 

With  the  use  of  explosive  sound  sources  with  a  broad  frequency  band, 

7 

Pekeris  found  that  the  effects  of  shear  waves  on  the  propagation  of  sound  in 

the  water  increases  with  decreasing  depth  of  the  water  and  with  decreasing 

frequency  of  the  propagated  sound.  Since  the  attenuation  coefficient  of  sound 

2 

is  directly  proportional  to  its  frequency,  this  and  the  low-frequency 

absorption  caused  by  the  shear  waves  makes  the  ocean  wave  guide  a  band-pass 

8 

filter  of  the  sound  emitted  by  the  source. 

If  the  acoustic  properties  of  this  ocean  wave  guide  are  well  known,  then 

9  10 

it  is  possible  to  calculate  the  optimum  frequency  of  sound  propagation.’  It 
has  been  observed  that  sound  with  frequency  below  the  optimum  frequency  highly 
depends  on  the  acoustic  properties  of  the  bottom  sediments  because  of  the 
relatively  small  absorptive  effects  of  low  frequency  sound,  while  the 
transmission  of  sound  with  frequency  greater  than  the  optimum  frequency 
depends  mostly  on  the  acoustic  properties  of  the  water  column!' 

Since  the  optimum  frequency  of  most  ocean  environments  is  smaller  than 
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one  kilohertz,  an  underwater  acoustic  model  for  the  propagation  of  sound  in 

this  frequency  band  is  needed.  This  model  must  contain  the  effects  of  the 

acoustic  properties  of  the  water  column  and  the  seismo-acoustic  properties  of 

1  2 

the  bottom  sediments. 

The  ray  theory  is  highly  satisfactory  to  predict  and  explain  some 
electromagnetic  phenomena,  and  it  is  very  useful  in  predicting  the 
transmission  loss  of  high-frequency  sound.  However,  this  asymptotic 
approximation  cannot  explain  wave-type  phenomena  such  as  diffraction  and 

13  14 

interference.  Pedersen  pointed  out  that  approximating  the  ocean’s  sound 
speed  profile  by  layers  of  constant  sound  speed  gradients  causes  erroneous 
transmission  loss  computations  where  acoustic  interference  occurs.  However, 
his  transmission  loss  calculations  are  made  using  ray  theory  which  is  the 
greater  cause  of  errors. 

The  normal  mode  theory  represents  the  exact  solution  to  the  perfect  wave 

guide  and  it  is  used  mostly  in  low-frequency  underwater  sound  predictions!^ 

This  model  could  also  be  used  to  predict  the  behavior  of  the  transmitted 

high-frequency  sound  with  the  expense  of  more  computation  time  and  memory. 

With  some  assumptions  and  approximations,  the  normal  mode  theory  can  be 

expanded  to  include  the  effects  of  layers  with  linear  acoustic  properties, 

shear  waves  from  elastic  layers,  and  range  dependence  of  the  ocean  wave  guide. 

Various  underwater  acoustic  models  have  been  developed  which  treat  some 

of  these  properties.  Each  model  has  its  virtues  and  limitations.  For  high- 

frequency  sound  propagation  the  ray  theory  can  handle  all  the  mentioned 

properties  but  we  are  interested  in  the  low  frequency  region. 

The  parabolic  equation  (PE)  method  was  originated  by  Frederick  Tappert'* 

in  the  late  70’ s  where  he  approximated  the  wave  equation  in  order  to  obtain  a 

parabolic  differential  equation  which  has  the  mathematical  property  of  a 

closed  form  solution  which  can  be  solved  by  calculating  the  transmission  loss 

while  incrementing  in  range  without  iterations.  This  approximation  involves 

neglecting  the  incoming  solution  of  the  wave  equation  and  the  use  of  the 

large-range  asymptotic  approximation  of  the  Hankel  functions.  This  method  was 

1  7 

adopted  by  Jensen  and  Kuperman  to  study  the  propagation  of  sound  in  an 
up-slope  wedge-shaped  wave  guide  with  a  fluid-type  bottom.  Their  results 
showed  that  as  the  depth  of  the  water  gets  shallower,  the  trapped  modes  reach 
their  cutoff  frequency  and  become  part  of  the  continuous  spectrum.  Later,  Ding 
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18  19 

Lee  ’  used  an  implicit  finite-difference  model  to  solve  the  parabolic 

20 

equation.  McDaniel  made  comparisons  between  normal  mode  and  the  parabolic 
approximation  concluding  that  the  latter  displays  errors  in  the  computation  of 
the  phase  and  group  velocities.  The  limitations  of  the  parabolic  approximation 
are  that  it  applies  only  to  sound  propagation  at  small  angles  with  respect  to 
the  horizontal,  it  does  not  handle  the  horizontally  reflected  waves,  and  it  is 
incapable  of  handling  the  shear  waves  created  in  the  elastic  bottom 

2122  23-27 

layers.  ’  Many  acousticians  have  included  shear  waves  in  their  PE 

models,  but  their  theory  is  based  on  approximations  and  lacks  the  rigorous 
derivation  of  elastic  effects. 

2  8 

The  fast  field  method  developed  by  Schmidt  is  made  to  include  the 
elasticity  of  the  bottom  layers,  but  it  may  only  be  used  in  range-independent 
environments . 

29  30 

Porter  and  Reiss  ’  have  also  included  the  shear  waves  by  incorporating 

the  impedance  condition  as  a  boundary  condition  to  the  normal  mode 

computations.  However,  this  is  also  a  range-independent  model  which  uses  the 

plane-wave  approximation  to  calculate  the  acoustic  impedances. 

The  coupled  normal  mode  method  was  implicitly  originated  by  Allan 
31  3  2 

Pierce  ’  in  the  mid  60s  as  an  adiabatic-mode  theory,  with  eigenray 

calculations  to  estimate  the  coupling  coefficients,  to  simplify  the  solution 

to  the  propagation  in  a  range-dependent  environment.  In  this  adiabatic  mode 

theory,  he  assumes  an  isovelocity  wave  guide  and  a  weak  coupling  between  the 

natural  modes  of  the  wedge-like  wave  guide.  He  found  that  compressional  waves 

33-35 

refract  into  the  basement  until  they  get  completely  attenuated.  McDaniel 
used  the  coupled  wave  equations  to  calculate  the  energy  transferred  between 
normal  modes  as  a  result  from  bottom  scattering  of  the  ocean,  and  has  shown 
that  randomly  rough  sea  bed  layering  can  increase  the  transmission  loss 
depending  upon  the  degree  of  penetration  of  the  acoustic  field  into  the 

3  6  —  38 

sediment.  Evans  modeled  the  axisymmetric  range-dependent  medium  as  N 

range- independent  segments  with  a  pressure-release  false  bottom  suitably  deep 
to  convert  the  continuous  spectrum  into  a  discrete  form.  The  eigenvalues  and 
eigenfunctions  of  each  range- independent  segment  was  solved  by  taking  into 
account  only  the  absorption  of  the  basement  layer  to  avoid  the  reflected 
energy  from  the  pressure-release  false  bottom.  The  group  made  up  of  by  Graves, 
Chwieroth,  Miller,  Nagl,  Uberall,  and  Zarur^^  have  used  a  similar  method 
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for  solving  the  set  of  coupled  range  equations,  but  they  included  water  layers 
with  linear  pressure  wave  number  squared.  This  gives  a  betters  approximation 
of  the  sound  speed  profile  and  the  solution  for  each  layer  is  given  by  Airy 
functions.  Comparisons  of  the  PE  model  and  the  coupled  normal  mode  model  have 

44  45 

been  made  for  some  up-slope  range-dependent  wave  guides  ’  which  gives  the 
best  agreement  if  the  environment  slowly  varies  with  range.  However,  it  has 
not  been  possible  to  define  how  slow  the  range  dependent  wave  guide  must  be. 

Even  though  the  elastic  effects  were  not  incorporated  in  the  coupled 
normal  mode  method,  there  is  no  reason  why  this  method  cannot  be  applied  to 
the  elastic  wave  equation.  The  objective  of  this  report  is  to  expand  the 
normal  mode  theory  in  order  to  predict  the  transmission  loss  of  low-frequency 
underwater  sound  with  the  effects  of  a  depth-dependent  sound  speed  in  the 
water  column  and  the  elastic  effects  of  the  solid  sediments  of  the  bottom. 
Since  some  range  dependence  of  the  acoustic  properties  and  the  boundaries 
between  the  layers  has  been  experimentally  observed,  this  property  will  be 
included  using  an  adiabatic  approach  assuming  slowly  range-dependent  acoustic 
properties.  Another  feature  of  the  model  to  be  presented  here  is  that  the 
unrealistic  false  boundary  used  by  Evans  and  Miller  to  convert  the  continuous 
spectrum  into  a  discrete  form  has  been  removed  because  the  absorptive  effects 
of  the  elastic  layers  causes  the  radiating  spectrum  to  be  discrete. 

The  derivation  of  the  wave  equation  and  theory  of  normal  modes  in  a 
liquid  wave  guide  is  developed  in  the  second  chapter.  Since  the  attenuation 
coefficient  of  low-frequency  sound  in  water  is  negligible  compared  to  the  loss 
in  the  bottom  elastic  sediments,  this  property  will  not  be  included.  However, 
each  layer  of  the  horizontally  stratified  water  column  will  have  a  constant 
density  and  sound  speed  gradient. 

In  the  next  chapter,  the  general  elastic  wave  equation  will  be  derived 
and  solved  by  dividing  the  ocean  floor  into  layers  of  constant  acoustic 
properties  and  separating  the  solution  into  a  divergent  (pressure)  and  a 
rotational  (shear)  term.  Each  solid  sediment  will  be  described  by  a  layer 
thickness,  a  density,  a  compressional  sound  speed,  a  compressional  attenuation 
coefficient,  a  shear  sound  speed,  and  a  shear  attenuation  coefficient. 

The  liquid-solid  boundary  conditions  are  derived  in  the  next  chapter  in 
order  to  match  the  solutions  for  each  adjacent  layer.  The  ocean  wave  guide  has 
a  set  of  trapped  modes  which  are  evanescent  in  the  bottom  and  a  set  of 
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radiating  modes  which  are  the  major  contributor  to  the  sound  propagation  into 
the  bottom  of  the  ocean  and  they  become  damped  by  absorption. 

A  complex  characteristic  equation  is  derived  where  the  complex 
eigenvalues  can  be  found  by  searching  for  the  complex  zeros.  The  most 
challenging  part  of  this  normal  mode  method  is  to  find  the  best  method  of 
determining  these  complex  roots  of  the  characteristic  equation,  since  "There 
are  no  good,  general  methods  for  solving  systems  of  more  than  one  nonlinear 
equation" The  Newton-Raphson  Method  for  Nonlinear  Systems  of  Equations^^  and 
the  Muller  Method  with  Deflation  used  in  the  subroutine  "DZANLY"  from  the  IMSL 

4  7 

Library  were  used  to  converge  into  the  complex  zeros,  but  severe  divergence 
has  been  experienced  for  some  of  the  roots.  Therefore,  the  Levenberg-Marquardt 
minimization  algorithm  for  the  magnitude  of  the  complex  determinant  has  been 
adopted  for  the  uniform  convergence  to  the  nearest  minimaf^ 

The  adiabatic  approximation  for  slowly  range-dependent  ocean  wave  guides 
is  finally  used  to  obtain  the  range-dependent  sound  propagation  model  with  the 
effects  of  the  elastic  and  absorptive  properties  of  the  ocean  floor.  With  such 
a  model,  it  will  be  possible  to  better  understand  the  various  propagation 
phenomena  related  to  low-frequency  sound  and  to  predict  the  transmission  loss 
of  sound  in  realistic  ocean  conditions. 
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CHAPTER  2 

DERIVATION  OF  THE  WAVE  EQUATION  FOR  FLUID  LAYERS 


The  wave  equation  is  a  mathematical  description  of  the  reaction  of  the 
media  due  to  a  disturbance  from  an  external  force  caused  by  a  source  or 
sources.  The  media  can  be  in  the  state  of  gas,  liquid,  or  a  solid,  and  the 
source  may  be  electromagnetic  or  mechanical.  In  this  chapter,  the  mechanical 
(acoustic)  propagation  of  the  disturbance  in  a  fluid  medium  is  treated.  The 
fluid  medium  is  the  ocean  environment  modeled  as  a  horizontally  stratified 
acoustic  wave  guide  where  the  surface  is  treated  as  a  pressure-release,  or 
soft,  or  resilient,  boundary.  The  bottom  is  modeled  as  elastic  layers  in  the 
next  chapter. 

The  disturbance  created  by  the  acoustic  source  may  be  expressed  as  a 
change  in  the  total  pressure,  relative  to  the  undisturbed  pressure,  as  a 
function  of  the  density  fluctuation  created  by  this  external  force.  If  the 
density  fluctuation  is  much  smaller  than  the  undisturbed  density  of  the 
environment,  then  the  total  disturbed  pressure  may  be  expanded  in  the 
following  Taylor  series: 


P(p)  =  P^+ 


dP 


I  5P  J 


(p-p  )  +  - 
0  2 


5^P 


ap‘ 


(p-p. 


(1) 


where  the  partial  derivatives  are  constants  determined  for  the  adiabatic 
compression  and  expansion  of  the  fluid  about  its  equilibrium  density  p^,  the 
equilibrium  pressure  is  P^,  and  the  instantaneous  total  density  is  p. 

If  the  magnitude  of  the  condensation  is  much  smaller  than  unity,  i.e.. 


s  H  (p  -  p^)/p^  =  pyp^ 


(2) 
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then  the  first  two  terms  in  the  Taylor  expansion  are  of  greatest  contribution 
and  an  acoustic  pressure  caused  by  the  disturbance  may  be  defined  as 


p  =  P(p)  -  s 


(3) 


4  9 

where  by  thermodynamic  arguments  it  is  found  that,  in  an  adiabatic  media, 
the  sound  speed  is  given  by 


2 


C 


'  dP  ' 


(4) 


and  the  adiabatic  bulk  modulus  is  given  by 


therefore,  the  acoustic  pressure  is  simplified  to 


2 

p  =  C 


(5) 


(6) 


which  is  called  the  state  equation. 

An  equation  for  the  motion  of  the  particles  in  the  fluid  is  also 
necessary  for  the  proper  environmental  description.  Consider  an  infinitesimal 
cubic  volume  in  the  medium  where  the  disturbance  is  taking  place  as  shown  in 
Figure  1(a)  for  the  one-dimensional  derivation  in  cartesian  coordinates. 
Equating  forces  in  a  continuous  medium  gives 


F  +  A  [P(x)  -  P(x+dx)]  =  ^(p  A  dx 

external  dt  dt 


(7) 
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where  the  external  force  is  the  disturbance  created  by  the  sound  source  and 
can  be  written  in  terms  of  a  "force  density"  with  the  expression 


F  s  X  A  dx  (8) 

external  e 


and  using  the  definition  of  a  derivative 


^  ^  P(x-t-dx)-P(x)  I 
5x  dx 


(9) 


Equation  (7)  becomes 


X 


5P 

9x 


V  |-(V  p)  +  |-{V  p) 

X  ox  X  ot  X 


which  in  three  dimensions  is  given  by 


(10) 


X  -  7P  =  V  (y-pV)  +  IrCpV)  (11) 

e  ot 

where  the  density  is  a  function  of  space  and  time,  and  the  equation  is  in  a 
non-linear  form.  The  total  pressure  is  P  and  the  total  instantaneous  particle 
velocity  is  V.  Dividing  the  instantaneous  density,  pressure,  and  particle 
velocity  into  an  undisturbed  part  and  an  acoustic  part.  Equation  (11) 
simplifies  to  the  linear  form 


5p  *  (12) 

where  the  acoustic  pressure  is  p  and  the  particle  velocity  is  v. 

Since  the  fluid  of  interest  is  continuous  throughout  the  infinitesimal 
volume.  Figure  1(b)  will  be  helpful  in  deriving  a  continuity  equation  under 
the  argument  that  the  mass  moving  into  the  volume,  p(x)  A  V^(x)  dt,  must  be 
the  same  as  the  mass  coming  out,  p(x+dx)  A  V  (x+dx)  dt.  There  may  be  a  change 

X 
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in  mass  inside  the  volume  due  to  the  compressibility  of  the  fluid,  ^  A  dx  dt, 

at 

and  there  may  exist  a  source  of  mass  inside  the  volume  represented  by  Q  A  dx 
dt.  Taking  the  definition  of  the  derivative  in  Equation  (9)  gives 


-  ^(pV)  A  dx  dt  =  ^  A  dx  dt  +  Q  A  dx  dt 
dx  St 

which  is  rewritten  in  three  dimensions  as 


(13) 


(pV)  +  ^  +  Q  =  0 


(14) 


or  in  a  linearized  form  as 


where  Q  =  0  when  the  source  is  external. 

Substituting  Equation  (6)  into  Equation  (15)  for  the  acoustic  density, 
taking  its  partial  derivative  in  time,  and  dropping  the  subscript  0  of  the 
undisturbed  density  of  the  medium,  gives 


P 


l^(V.v) 


-2 

+  C 


and  taking  the  divergence  of  Equation  (12)  yields 


(16) 


p7- (1/p  Vp) 


^  P^(V-v)  = 


p  V- (X/p) 


(17) 


which  subtracted  from  Equation  (16)  provides  the  inhomogeneous  wave  equation 


p7-(l/p  Vp)  +  kS  =  P  V-(l/p  W) 


(18) 
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where  the  external  force  has  been  written  in  terms  of  an  external  potential 
energy,  time  harmonic  behavior  has  been  assumed  where  k=u/c,  and  the 
undisturbed  density  of  the  fluid  is  taken  as  space  dependent.  This  equation 
can  also  be  written  as 

y^p  +  -  p'\vp)-ivp)  =  -  p‘^(yp)-(yy)  (19) 

which  is  simplified  under  the  change  of  variables 


p  =  n  (20) 

and 

U  =  Vp  V  (21) 

to  obtain 

y^n  +  (k^+K^)n  =  (22) 

where 

(23) 

If  the  density  is  taken  as  a  linear  function  of  deptn,  then  Equation  (23) 
simplifies.  However,  the  inhomogeneous  equation  to  solve  also  has  a  depth 
dependent  wave  number  to  worry  about  due  to  the  depth  dependence  of  the  sound 
speed.  The  changes  in  density  with  depth  only  occurs  in  the  surface  of  the 
Arctic  Ocean  where  water  has  been  solidified  and  at  regions  where  the  water 
from  rivers  merge  into  the  salty  ocean  waters.  It  is  concluded  that,  for 
simplicity,  the  ocean  environment  can  be  divided  into  horizontal  layers  with 
constant  density.  It  is  understood  that  the  bottom’s  solid  sediments  may  have 
layers  of  large  density  gradients?^  In  this  case.  Equation  (22)  must  be 
solved.  However,  it  is  also  possible  to  divide  this  layer  into  smaller 
constant-density  layers.  Then  Equation  (18)  becomes 

y^p  +  k^p  =  V^U  (24) 
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which,  in  an  unbounded  medium,  has  a  general  homogeneous  solution  consisting 
of  an  outgoing  and  an  incoming  wave,  and  an  inhomogeneous  solution  caused  by 
the  external  force.  Since  a  sound  source  in  a  fluid  can  only  produce  a  scalar 
potential  (no  shear  waves),  the  curl  of  Equation  (12)  gives  the  property 


7xv  =  constant  s  vorticity  (25) 

the  vorticity  in  the  medium  does  not  change.  Therefore,  if  initially  there 
has  been  no  rotational  component  of  the  particle  velocity  then  the  vorticity 
will  always  be  null  and  this  particle  velocity  can  be  written  in  terms  of  a 
velocity  potential 


V  =  V(p 


(26) 


which  substituted  back  into  Equation  (12)  gives 


P  =  -P 


d<p 

at 


(27) 


and  assuming  harmonic  time  dependence  yields 


p  =  -Uj}p  <f> 


(28) 


or 


V  = 


-■cVp 


up 


(29) 


Substitution  of  Equation  (28)  into  Equation  (24)  provides  the  inhomogeneous 
Helmholtz  equation 


n2  ,2  <.  -2,, 

V  tp  +  k  w  =  -  7  U 

^  (j  p 


(30) 


which  must  be  solved  for  the  velocity  potential.  If  the  medium  is  bounded. 
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the  solution  must  satisfy  the  appropriate  boundary  conditions. 

The  conservation  of  energy  is  obtained  by  the  scalar  product  of  Equation 
(12)  with  the  particle  velocity  and  substituting  the  continuity  equation, 
Equation  (15),  providing  the  law  of  conservation  of  energy: 


where 


dS  p.  ■* 


(31) 

(32) 


is  the  acoustic  energy  density,  and 


1  =  p  V 


(33) 


is  the  acoustic  energy  flux  or  acoustic  intensity  in  SI  units  of  Watts  per 
square  meter.  Integrating  Equation  (31)  throughout  a  volume  in  the  fluid 
medium  provides  the  acoustic  power 


e  dV  =  -  i  I-n  dS  (34) 

S 

in  terms  of  a  closed  surface  integral  around  the  volume  where  all  the  energy 
is  contained. 

To  obtain  these  important  measurable  quantities,  it  is  necessary  and 
sufficient  to  solve  Equation  (30)  for  the  velocity  potential.  To  solve  this 
equation  by  separation  of  variables,  it  has  been  proven  that  the  sound  speed 
must  be  a  function  of  only  one  variable?^  This  variable  is  taken  to  be  the 
vertical  direction  since  temperature  and  the  total  pressure  of  the  ocean 
highly  depends  on  depth.  Based  on  experimental  measurements,  numerical  fits 
have  been  made  to  determine  the  sound  speed  and  its  attenuation  coefficient  in 
the  sea. 

A  simplified  version  of  Wilson’s  formula  for  the  sound  speed  as  a 
function  of  temperature,  salinity,  and  depth  is  given  by 
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c(z)  =  1492.9  +  3(T-10)  -  6xl0‘^(T-10)^  -  4xlO"^(T-18)^  +  1.2(5-35)  - 

10'^(T-18) (S-35)  +  z/61  (35) 

where  the  temperature  T  is  in  Celsius,  the  salinity  S  is  in  parts  per 
thousand,  and  the  depth  z  is  in  meters.  The  formula  is  accurate  to  0.1  m/s 
for  a  temperature  less  than  20°C  and  for  depths  less  than  8.0  kilometers^  A 
sound  speed  profile,  with  its  salinity  and  temperature  profiles,  is  displayed 
in  Figure  2  as  a  function  of  depth  in  the  Arctic  Ocean.  Note  the  variability 
of  the  data  which  is  caused  by  currents.  Sound  speeds  in  the  oceans  may  vary 
anywhere  between  1400  m/s  and  1600  m/s. 

The  attenuation  coefficient  has  been  fitted  by  Thorp  as  a  function  of 
frequency,  temperature,  and  pressure.  A  simplified  version  of  the  formula  is 

a(f)  =  0.003  +  0.1f^/(l+f^)  +  40fV(4100+f^)  +  0.000275f^  (36) 

where  f  is  the  frequency  in  kilohertz,  a  is  the  attenuation  coefficient  in 
decibels  per  kiloyard,  and  the  relationship  holds  for  frequencies  greater  than 
ten  hertz  but  lower  than  one  megahertz.  This  attenuation  coefficient  is 
included  in  the  environmental  properties  as  the  imaginary  part  of  the  wave 
number.  However,  it  must  be  converted  to  units  of  nepers  per  meter.  The 
conversion  is  given  by 

a(nepers/meter )  =  a(dB/kyd)  /  [20  log(e)]  /  914.4.  (37) 

Note  that  the  attenuation  coefficient  is  directly  proportional  to  the 
frequency,  therefore  low  frequency  sound  travels  farther  than  its  counterpart. 
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CHAPTER  3 

THE  GENERAL  ELASTIC  WAVE  EQUATION 

Consider  two  infinitesimal  volume  elements,  P  and  Q,  at  a  distance  Ar 
from  each  other  in  an  elastic  medium.  An  external  disturbance  moves  these 
elements  to  the  position  P’  and  Q’ ,  respectively. 

Strain,  a  dimensionless  quantity,  is  defined  as  the  change  in  position  of 
a  point,  say  Q,  with  respect  to  a  reference  neighboring  point,  say  P,  divided 
by  the  distance  between  these  points,  i.e.  /Ar/.  Therefore,  expand  each 
component  of  the  ^  displacement  in  a  Taylor  series  relative  to  the  ^ 
displacement  as  follows: 

5^  ac  ac 

<:  =  ?  +  Ax  +  Ay  +  -5-^  Az  +  ...  (38) 

X  ^x  ax  ay  dz 

and  combine  them  to  obtain  the  expansion, 

^  ^  +  (Ar  •  ^  .  (39) 

where  the  higher  order  terms  are  nonlinear  and  may  be  neglected  if  we  assume 
the  maximum  amplitude  of  the  displacement  vectors,  C  and  <!;,  to  be  much  smaller 
than  the  wavelength  of  the  disturbance  (Hooke’s  law).  This  assumption  is 
valid  at  low  frequencies  and  at  distances  greater  than  a  wavelength  from  the 
external  source. 

Adding  and  subtracting  the  terms  Az  and  ■  Ay  in  Equation  (38) 

gives, 

<  =  C  +  Az  -  n  Ay  +  S  Ax  +  S  Ay  +  S  Az  +  . . .  (40) 

X  X  y  z  XX  xy  xz 

with 
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and 


a  = 


1 

2 


^  xt 


(41) 


(42) 


where  t  stands  for  the  transpose  of  the  tensor  matrix  in  the  parenthesis  and 
where  <*  identifies  a  tensor  of  second  rank.  Although  there  are  no  units  for 
strain,  engineers  sometimes  use  implied  multiples,  e.g.  microinches/inch,  a 
number  which  is  strain  x  10^,  or  percent  strain,  a  number  which  is  strain  x 
10  .  Where  strain  is  displacement  per  unit  of  length  in  the  direction  of 
displacement,  it  is  referred  to  as  a  normal  strain.  Where  the  strain  is 
displacement  per  unit  of  length  in  a  direction  perpendicular  to  the  direction 
of  displacement,  it  is  referred  to  as  a  shear  strain.  In  a  three-dimensional 
strain  field  the  strain  can  be  resolved  into  an  isotropic  component  (or 
volumetric  strain),  representing  change  in  volume  (and  density)  of  an  element, 
and  a  deviatoric  component,  representing  change  in  the  shape  of  an  element. 

In  vector  form  we  get. 


^=^+^xAr+§  .  Ar+  ...  (43) 

where  ^  is  expressed  in  terms  of  the  rotational  and  divergent  behavior 
relative  to 

In  an  elastic  medium,  the  linearized  Euler  equation  becomes, 

+  p  =  ?  (44) 

where  P  is  the  stress  tensor,  v  =  d^/dt  is  the  particle  velocity,  p  is  the 
undisturbed  density  of  the  medium,  and  t  is  the  external  force  that  disturbs 
the  medium  creating  the  acoustic  field.  Stress  is  force  per  unit  area  and  may 
be  either  normal  stress  produced  by  tensile,  compressive  forces  acting 
perpendicular  to  the  faces  of  cubic  elements,  or  shear  stresses,  produced  by 
tangential  forces  acting  parallel  to  the  surfar*»s  of  cubic  elements.  In  a 
three-dimensional  stress  field,  the  stress  sy.  -m  can  also  be  resolved  into  an 
isotropic  component  (or  bulk  stress)  and  a  deviatoric  component.  By 
definition,  this  elastic  medium  has  a  rest  state  in  which  stresses  and  strains 
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vanish.  Linear  elasticity  means  that  stress  is  directly  proportional  to 
strain.  If  the  medium  is  also  isotropic  then  the  relationship  between  the 
stress  and  strain  tensors  in  terms  of  Lame  constants  is  given  by  the 
constitutive  relation 


-?  =  Xef  +  2/i§  (45) 

where  0  =  ^  ^  is  the  divergent  component  of  the  displacement  vector,  ?  is 

the  unit  matrix,  and  the  Lame  constants  are  given  by 


and 


X  =  p(c^  -  2  b^) 


(46) 


p  =  p  b 


(47) 


where  c  is  the  compressional  sound  speed  of  the  medium  and  b  is  its  shear 
speed.  Equations  (45)  through  (47)  are  sufficient  to  characterize  the  linear, 
isotropic,  elastic  medium.  The  absorptive  nature  of  the  medium  is  modeled  by 
making  both  sound  speeds  complex  where  the  imaginary  part  of  the  sound  speeds 
is  related  to  their  respective  attenuation  coefficients.  Since  the  speeds  are 
considered  space  dependent,  the  substitution  of  the  stresses  into  Euler’s 
equation.  Equation  (44),  gives 


and  substitution  of  Equations 


^  •  (X  0  ?)  +  2  ^  •  (p  §),  (48) 

(41)  and  (42)  provides  the  elastic  wave  equation 


p  ^  +  (X  +  p)  ^0  +  0  ^X  +  p  +  2  •  §  (49) 

5t^ 

which  is  a  set  of  three  coupled  inhomogeneous  differential  equations  where 
each  component  of  the  solution  ^  depends  on  the  others. 

To  uncouple  this  generalized  elastic  wave  equation,  it  is  sufficient  to 
eliminate  the  last  term  in  the  right-hand-side.  This  is  done  by  assuming  a 
quasi  space-independent  shear  modulus,  p.  From  Equation  (47)  this  means  layers 


16 


NSWC  TR  89-170 


of  constant  shear  speed  and  density,  but  note  that  no  restriction  is  imposed 
on  the  compressional  speed. 

In  the  case  of  solid  (elastic)  layers  there  is  no  external  force  present, 
?  =  5,  since  the  source  is  assumed  to  be  in  the  water  column.  Furthermore, 
layers  of  constant  shear  modulus  are  assumed.  Under  these  conditions  Equation 
(49)  simplifies  to 


p  I"  =  (X  +  fi)  +  6 


(50) 


Since  the  displacement  vector  has  been  expanded  into  its  rotational  and 
divergent  part  in  Equation  (43),  the  same  expansion  may  be  applied  to  the  time 
derivative  of  the  displacement  vector  now  defined  as  the  particle  velocity, 
i.e.  , 


V  =  ^  y  t 


(51) 


where  the  first  term  in  the  right-hand-side  corresponds  to  the  compressional 
waves  and  the  second  term  is  the  shear  wave  contribution  to  the  acoustic 
field.  Substituting  Equation  (51)  into  Equation  (50),  and  taking  into 
consideration  the  assumption  of  isodensity  layers  with  constant  shear  speed, 
gives 


^  p  -  (A  +  2m)  j  =  ^  X  M  -  P  •  (52) 

The  divergence  of  this  vector  equation  is  satisfied  if 

(V^  +  k^)  (p  =  0,  (53a) 

which  is  the  compressional  wave  equation,  and  the  rotational  part  is  satisfied 
with  the  shear  wave  equations 


(V^  +  K^)  X  =  0, 


(53b) 
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where  k  =  w  /  c  is  the  compress ional  wave  number,  k  =  w  /  b  is  the  shear  wave 
number,  and  harmonic  time  dependence  is  assumed. 

For  simplicity,  the  compressional  and  shear  attenuation  coefficients  are 
included  as  the  imaginary  part  of  both  wave  numbers  instead  of  the  sound 
speeds.  In  this  case,  the  complex  compressional  wave  number  is 

k  =  u  /  c  *  i  a,  (54a) 

where  a  is  the  compressional  attenuation  coefficient  and  the  complex  shear 
wave  number  is 


K  =  u  /  h  +  i  8,  (54b) 

where  8  is  the  constant  shear  attenuation  coefficient. 

For  the  case  of  axially  symmetric  propagation  in  cylindrical  coordinates 
we  write 


V  =  V  (r,z)  r  +  V  (r,z)  z  (55) 

r  2 

which  will  be  equated  to  Equation  (51)  to  obtain  the  components  of  the 
particle  velocity  in  terms  of  the  potentials.  However,  to  obtain  range 
independent  boundary  conditions  the  vector  potential  A  must  be  written  as  the 
curl  of  another  vector  potential,  i.e.. 


t  =  ^  X  '$  (56) 

which  gives 

V  =  +  ^  (^  •  i^)  -  (57) 

The  physical  meaning  of  this  new  vector  potential  is  obvious  when 
calculating  the  stress  matrix  for  the  boundary  conditions.  The  component  {p 

5  2  ^ 

corresponds  to  the  SH-waves  in  seismology,  0  vanishes  in  an  axially 

0 

symmetric  medium,  and  corresponds  to  the  SV-waves.  However,  the  stress 
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tensor  of  the  SH-wave  has  only  off-diagonal  elements,  while  the  stress  tensor 
in  the  liquid  layers  has  null  off-diagonal  elements.  Therefore,  a  source  in 
the  liquid  layers  is  incapable  of  exciting  SH-waves  in  the  solid  layers.  In 
consequence,  we  have  tp  =  0  and  the  particle  velocity  is  given  by. 


■  a  a 

3  , 1  ■  ^ 

r  3 

1  a 

r  a  .  11 

V  = 

dr  ^  *  dz 

r  ar 

^ - 

^1 

N 

Z 

The  shear  wave  equations  to  be  satisfied  are  now  reduced  to  the  scalar  wave 
equation. 


(7^  +  0  =  0. 


Substituting  Equations  (59)  into  Equation  (58),  equating  this  result  with 
Equation  (55),  and  assuming  axially  symmetric  cylindrical  coordinates  provides 


a 

V  =  — 

2  8z 


a  r  a^ 


+  (C  \  [(I 


(60a) 


(60b) 


for  the  components  of  the  particle  velocity. 
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CHAPTER  4 

SOLUTION  OF  THE  WAVE  EQUATION 


The  solution  of  the  inhomogeneous  compressional  Helmholtz  equation  for 
the  liquid  layers  can  be  written  as  the  sum  of  the  homogeneous  solution  and 
the  particular  (transient)  solution.  The  generalized  homogeneous  solution  can 
be  used  as  the  solution  of  the  homogeneous  compressional  wave  equation  in  the 
solid  layers,  Equation  (53a).  The  solution  of  the  homogeneous  shear  wave 
equation,  Equation  (59),  however,  is  different  since  the  shear  speed  and 
compressional  speed  of  a  sediment  are  usually  unequal. 

For  simplicity  and  without  loss  of  generality,  we  may  solve  the 
inhomogeneous  wave  equation  in  the  water  column  for  the  case  of  a  point  source 
of  unit  source  strength  which,  in  the  cylindrical  coordinate  system,  becomes 

2 

:j  +  +  k^(z)j  ^(r,z)  =  -  5(r)  6(z  -  z^)  (61) 

where  the  source  isatr=0,  z=z,  and  the  wave  number  is  taken  to  be  depth 

O 

dependent  only. 

A  lengthy  but  elegant  way  to  solve  the  wave  equation  is  by  the  use  of  the 
Fourier-Bessel  transformation. 


i  ^ 

r  dr 


<p{r,z) 


u(k,z)  k  dk. 


and  its  inverse  form. 


u(k, z) 


(fir.z)  Jp(ki') 


r  dr. 


(62a) 


(62b) 
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where  the  zeroth  order  Bessel  functions  satisfy  the  closure  relation, 


6(r  -  r’ ) 


J  (kr)  J  (kr'  ) 
0  0 


k  dk. 


(63) 


and  the  Bessel  equation, 


k^r^  J"(kr)  +  kr  J’ (kr)  +  k^r^  J  (kr)  =  0. 
0  0  0 

Substitution  of  the  above  equations  and  the  relation 


(64) 


J  (kr)  =  k  J’ (kr)  +  k^r  J"{kr)  =  -k^r  J  (kr)  (65) 

0  0  0  0 

into  Equation  (61)  converts  the  partial  differential  equation  to  the  ordinary 
differential  equation 


,2  ^  T  6(z  -  z  ) 

-  +  k^(z)  -  k  1  u(k,z)  =  -  - J — — 


(66) 


2 

where  k  is  the  continuous  eigenvalue  and  u(k,z)  is  the  continuous 
eigenfunction  of  the  inhomogeneous  equation.  In  the  case  of  a  discrete  wave 
number  spectrum  the  eigenequation  becomes 


r  ,2  1 

-  +  k^(z)  -  k^  u  (z) 

,2  n  n 

I-  dz  J 


6(z  -  Z  ) 
_ 0 

Zn 


(67) 


2 

where  n  =  1,  2,  3,  . . . ,  N  is  the  mode  index,  k  are  the  discrete  eigenvalues 

n 

and  u  (z)  are  the  discrete  eigenfunctions.  These  eigenequations  are  similar  to 
n 

the  Schrddinger  equation  in  quantum  mechanics.  A  sound  speed  profile  taken  in 


21 


NSWC  TR  89-170 


the  Arctic  is  shown  in  Figure  2.  The  sound  speed  as  a  function  of  depth  is 
2  2  2 

contained  in  k  (z)  =  u  /c  (z)  and  the  minimum  sound  speed  is  defined  as  c 

Bln 

which  represents  a  maximum  wave  number  in  the  eigenequation  k  ,  Define  the 

max 

2  2 

"equivalent  potential"  as  V(z)  2  k  -  k  (z)  and  the  "equivalent  eigenvalues" 
2  2 

as  E  2  k  -  k  ,  then  the  eigenequation  becomes 
n  max  n 


u"(z)  +  [E  -  V(z)]  u  (z)  =  -  6{z  -  z  )/2n  (68) 

n  n  n  0 

which  is  equivalent  to  the  time-independent  Shrddinger  wave  equation,  where 

the  double  prime  stands  for  the  second  derivative  with  respect  to  the 

argument,  the  minimum  potential  is  zero,  and  there  exists  a  maximum  potential 

which  represents  the  threshold  between  the  discrete  and  continuous  spectrum. 

This  maximum  potential  is  given  by  the  determination  of  the  sound  speed  in  the 

limit  as  z  — >  «.  The  "propagating  modes"  are  defined  by  the  mode  cutoff  limit 
2  2 

E  <  u  /c  and  they  may  be  "trapped  modes'  which  are  represented  by  a 
n  Bin 

discrete  eigenvalue,  or  they  may  be  "radiating  modes"  which  are  represented  by 
the  continuous  spectrum.  Note  th^t  the  potential  V(z)  and  the  eigenvalues  E 

n 

depend  on  the  frequency  of  the  source.  Therefore,  the  number  of  discrete 
eigenvalues  vary  with  the  frequency  of  the  sound  that  is  disturbing  the 
medium. 

Since  the  ocean  floor  is  not  rigid  nor  soft,  the  energy  spectrum  will 
contain  trapped  modes  which  are  evanescent  in  the  bottom  layers  with  higher 
sound  speed  and  radiating  modes  which  represent  the  energy  that  radiates  into 
the  bottom  sediments. 

Evans^^  and  Miller*^"^^  have  solved  the  purely  real  wave  equation 
which  gave  a  discrete  eigenvalue  spectrum  representing  the  trapped  modes  and  a 
continuous  spectrum  representing  the  radiating  modes.  The  absorptive  effects 
of  the  fluid-type  bottom  was  later  incorporated  as  the  imaginary  part  of  the 
eigenvalues  using  the  first  order  perturbation  approximation  assuming  small 
absorptive  effects.  The  continuous  spectrum  was  forced  to  be  discrete  by 
adding  a  deep  false  (pressure-release)  boundary.  Since  this  boundary  caused 
reflection  of  the  incident  sound,  a  very  large  attenuation  coefficient  must  be 
used  in  the  basement  layer,  therefore  making  the  first  order  perturbation 
approximation  an  invalid  method. 
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The  continuous  eigenvalue  spectrum  is  the  direct  consequence  of  solving 
the  unrealistic  purely  real  wave  equation.  If  the  absorption  is  incorporated 
in  the  wave  equation  as  the  imaginary  part  of  the  wave  number,  then  the 
radiating  modes  will  be  part  of  a  discrete  spectrum.  In  this  case,  the 
Helmholtz  equation  in  Equation  (67)  becomes 


r  ^ —  +  k^(z)-k^]  u(z)=0  (69) 

L  dz"  "  J  " 

which  is  the  homogeneous  Helmholtz  wave  equation. 

From  the  continuity  of  pressure  in  the  liquid  layers,  which  will  be 
discussed  in  the  chapter  on  the  boundary  conditions,  we  shall  set  the  function 
v^P  (z)  u  (z)  as  the  orthonormal  eigenfunctions, 

n 


I 


p(2) 


u  (z)  u  (z)  dz  =  5 

n  m  nm 


(70) 


which  also  satisfies  Sturm-Louville’ s  Theorem?'  The  closure  relation  is  given 
by 


5(z  -  z  )  =  p(z  ) 
0  0 


N 


n=l 


(71) 


and  the  eigenfunction  is  solved  by  the  method  of  separation  of  variables  for 
the  individual  modes  which  are  later  added  in  the  form 


u(k,z) 


n 


n  =  l 


(k) 


u  (z) 

n 


(72) 


where  a  (k)  is  to  be  determined.  The  inhomogeneous  term  is  taken  into  account 
n 

if  we  substitute  the  homogeneous  solution  in  the  inhomogeneous  equation.  The 
substitution  gives 
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f  ^ —  +  k^(z)  -  1  ^  a  (k)  u  (z)  =  -  — \  u  (z)  u  (z  )  (73) 

I  dz^  i  ^  "  2n  n  no 

n=l  n=l 


and  substituting  the  homogeneous  equation  provides 


p(z  ) 

2  ,  2 ,  0 


N 

na  (k)  (k‘  -  k‘)  -  -  ^ 

n  n  2ti  n  0 

n=  1 


u 

n 


u  (z  )  u  (z)  =0 


(74) 


and,  to  satisfy  the  equation,  the  terms  inside  the  brackets  are  set  to  zero 
leading  to  the  relation 


p(z)  u  (z  ) 

/■1  1  0  n  0 

a  (k)  =  ^  ; - - 

"  ^  k  -  k^ 


(75) 


which  substituted  into  Equation  (72)  gives 


u(k,z)  = 


p(z  )  r — .  u  (z  )  u  (z 


/  t  »  U  \  ^  J  U  \ 

0  \  n  0  n 

2"  Z_  -  k" 

n  =  l  n 


(76) 


and  this  substituted  into  Equation  (66)  gives  the  scalar  potential 


n=l 


foo  J  (kr)  k  dk 
0 _ 

k  -  k 

n 


(77) 


where  the  integral  of  this  equation  is  better  solved  by  contour  integration 
and  it  can  be  defined  as 
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‘oo  H^^kr)  +  (kr) 

— - - -  k  dk.  (78) 

1  2  ,2 

k  -  k 

n 

A  property  of  the  eigenvalues  of  the  problem  is  that  these  usually  have  a 
smaller  imaginary  component  compared  to  the  real  part.  Also,  both  parts  of  the 
eigenvalues  are  positive  because  outgoing  waves  from  the  source  are  used  to 
satisfy  causality.  To  solve  this  complex  integral,  consider  the  contour 
integration  in  the  first  quadrant  of  the  complex  k-plane  as  displayed  in 
Figure  3,  where 


I  (r)  =  1 

n  2 


I 

ni 


(r) 


H  (kr ) 
0 _ 

1,2 

r  -  ^  ~ 

C  n 

1  i 


(kr) 

0 


k  dk. 


21 


k^-  k^ 

n 


k  dk. 


1  =  1,2,3 


1  =  4,5,6 


(79) 


and,  by  this  definition,  the  integral  to  be  solved  is. 


I  (r)  =  I  (r)  +  I  (r)  1. 

ni  2  L  n4  J 

4  9 

By  Jordan’ s  lemma  we  have 


(80) 


I  (r)  =  I  (r)  =  0 

n2  n5 


(81) 


also  the  integrals  in  Equation  (79)  show  that 


I  (r)  +  I  (r)  =  0 

n3  n6 


(82) 
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which  means  that  we  may  write 


I  (r) 


1 

2 


(r) 


(83) 


where  only  I  (r)  and  I  (r)  contribute  to  the  sum. 
nl  n4 

Given  that  the  singularities,  k  =  k  ,  are  located  in  the  upper  contour  we 

n 

get  that 


(kr) 
0 _ 

k^-  k^ 


k  dk 


0 


(84) 


and 


h“’  (kr) 
0 

k^-  k^ 

n 


k  dk  =  n  i  H*''(k  r) 

0  n 


(85) 


by  calculus  of  residues.  Substitution  of  Equations  (84)  and  (85)  in  Equation 
(83)  and  this  one  into  Equation  (77)  gives 


V(r, z) 


(z  )  u  (z)  h”’  (k  r ) 

On  On 


n=l 


(86) 


where  the  eigenfunctions  u  (z),  and  eigenvalues  k  ,  satisfy  Equation  (69). 

n  n 

53 

Note  that  the  solution  has  been  written  using  a  separation  of  variables. 

To  solve  the  characteristic  equation  for  the  compressional  waves, 
Equation  (69),  the  sound  speed  profile  is  divided  into  layers  where  the 
squared  of  the  index  of  refraction  is  a  linear  function  of  depth,  i.e., 
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n  (z)  *  a  z  +  b 
J  J  J 


where  k^(z)  =  w  n^Cz),  Figure  4  gives  the  geometry  to  be  used  for  this 
mathematical  model,  and  the  subscript  j  is  the  layer  number.  To  determine  a^ 
and  bj,  let  the  sound  speed  at  the  top  of  the  layer  to  be  c  and  that  of  the 
bottom  to  be  Substituting  into  our  linear  equation  gives 


-  =  a  2  +  b  , 

c2  J  J  J 


(88a) 


—  =  a  z  +  b  , 

c2  J  J*1  J 

J*1 


(88b) 


which  solved  for  a  and  b  results  in 

J  J 


2  2 
C  -  C 
J 

r>  2  2 

Dec 
J  j  J*1 


(89a) 


,2  2 

Z  (c  -  C 
J  J  J^l 

r,  2  2 

Dec 
J  J  J-1 


(89b) 


where  D  =  z  -  z  is  the  thickness  of  the  layer.  The  index  of  refraction  is 
J  J+i  J 

^  2  2  2  2 
given  by  n  =  1/c  ,  and  if  c  s  (n  -  n  )/D  ,  then  we  have 
®  J  J  J  J  J 


where 


=  n^  +  cr  (z  -  z)]  s  k^  +  S 
L  J  J  J  J  J  J 


(z  -  z  ) 
J 


(90a) 


S  5  (k^  '  k^  )/D 
J  J  J*i  J 


(90b) 
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and  which  substituted  into  the  eigenequation  gives 


dz 


u  (z)  +  [k^+S(z-z)-k^lu  (z)  =  0 

nj  L  J  J  J  n  J  nj 


(91) 


Define 


<  (2)  ^ 
n  J 


-  S  [  k^  +  S  (z  -  z  )  -  k^  1 
i  L  J  J  J  "  J 


(92) 


and  square  its  derivative  to  obtain 


d^  _  g  2/3  d^ 

dz"  ■  J  ^ 

nJ 


(93) 


which  substituted  into  the  new  eigenequation  gives 


[  — -  -  C  .  ]  u  (C  )  =  0.  (94) 

L  dC  -I 

nJ 

The  solutions  of  this  differential  equation  are  the  Bessel  functions  of  order 
1/3,  or  more  commonly  known  as  the  Airy  functions,  i.e., 


u 

nJ 


n  J 


a 

nJ 


Ai(<  )  +  b  Bi(c  )• 

nJ  nJ  ^nj 


(95) 


The  behavior  of  the  real  part  of  these  functions  and  their  derivative  is  shown 
in  Figure  5  and  some  of  their  properties  are  given  in  Reference  54.  Now  that 
the  general  solutions  are  found,  we  must  match  the  solutions  at  each  boundary 
with  the  appropriate  boundary  conditions  in  order  to  find  the  unknown 
coefficients  and  eigenvalues  for  each  mode. 

The  attenuation  coefficient  in  the  water  at  low  frequencies  is  negligibly 
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small.  Therefore,  only  the  shear  and  compressional  attenuation  coefficients  in 
the  elastic  bottom  layers  are  taken  into  consideration,  and  only  the  complex 
eigenvalue  will  make  the  argument  of  the  Airy  functions  complex. 

Using,  once  more,  the  method  of  separation  of  variables  for  the  solution 
of  the  homogeneous  shear  Helmholtz  equation  we  obtain  the  general  form 


0(r , z) 


V  (z)  H  (r) 

n  n 


which  substituted  back  in  the  shear  wave  equation  leads  to 


(96) 


1  d  d  1 


H  r  drl  dr 

n 


Ti  Id  2,2 

H  = - V  +  K  =  k  =  constant 

n  V  2  n  n 

n  dz 


(97) 


which  gives  the  ordinary  differential  equations. 


% 

-  V  (z)  +  [k^  -  k^]  V  (z)  =  0, 

d^2  n  n  n 

and 


2d  ^  d  ,22),.,, 

r  -  +  r  +  k  r  H  (r)  =  0. 

,2  dr  n  I  n 

'•  dz  ' 


The  solution  of  the  first  equation  in  the  layer  is 


(98) 


(99) 


V  (z)  =  c  exp(r  z)  +  d  exp(-y  z)  (100) 

nj  nj  ^  nj  nj  ^  nj 


2  2  2  2 
where  y  s  k  -  k  and  the  solution  can  be  oscillatory  (y  <  0)  or 

nJ  J  ^  '’nJ 

exponential  (y  >0).  By  causality,  only  the  radially  outgoing  solution  of 
n  J 

the  second  equation  must  be  used.  The  solution  to  this  radial  diffeiential 
equation  is  the  zeroth  order  Hankel  function  of  the  first  kind,  H*^*(k  r). 

On 

Finally,  since  a  and  b  are  unknowns  to  be  evaluated  by  the  use  of  the 
nj  nj 

boundary  conditions,  we  are  free  to  choose  A  s  i/A  p(z  )  u  (z  )  which,  in  the 

n  0  n  0 
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.th  , 

J  layer, 


gives 


(r.z) 


(z 

/  nj  0 


)  V  (z)  ’  (k  r) . 

nJ  On 


(101) 


n=l 


This  expression  simplifies  the  boundary  conditions  for  the  evaluation  of  the 
unknown  eigenvalues  and  amplitude  of  the  eigenfunctions. 

Note  that  the  n^^  compressional  eigenfunction  and  the  n*"^  shear 
eigenfunction  are  represented  by  the  same  eigenvalue  k  .  The  common  eigenvalue 

n 

is  necessary  to  satisfy  the  general  elastic  wave  equation  and  to  obtain 
range- independent  boundary  conditions. 

Now  that  the  general  solutions  are  found,  these  solutions  will  be  matched 
at  each  boundary  with  the  appropriate  boundary  conditions  in  order  to  find  the 
unknown  coefficients  and  eigenvalues  for  each  mode. 
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CHAPTER  5 

THE  BOUNDARY  CONDITIONS 


The  ocean  Is  modeled  as  a  horizontally  stratified  wave  guide  with  layers 
of  elastic  properties  simulating  the  ocean  floor  and  liquid  layers  simulating 
the  water  column.  To  match  the  solutions  of  adjacent  layers,  boundary 
conditions  are  derived.  In  this  chapter,  the  boundary  conditions  for  all 
possible  interface  are  developed. 


BOUNDARY  BETWEEN  LIQUID  LAYERS 

There  are  several  ways  of  calculating  the  boundary  conditions.  The 
boundary  conditions  for  the  interface  between  two  fluid  layers  can  be  obtained 
when  an  infinitesimal  cylindrical  volume  is  modeled  across  this  boundary. 

There  are  two  boundary  conditions  to  be  satisfied  at  this  interface. 


Continuity  of  the  Normal  Particle  Velocity 

The  volume  integration  of  Equation  (15)  in  this  infinitesimal  cylinder 
provides  the  expression 


v-n  dA  =  -  ^  J  p  dA  Ax  (102) 

where  making  Ax  — >  0,  the  right  hand  side  of  the  equation  vanishes  and  the 
surface  integral  yields  the  boundary  condition 


V  •  n  I  =  V  •  n  I  (103) 

2  ' s  1  ‘s 
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which  means  that  the  normal  component  of  the  particle  velocity  is  continuous 
at  the  interface. 

This  boundary  condition  is  expressed  as  v  •  n  =  continuous,  and  in  the 
case  of  horizontally  stratified  layers  we  may  write  n  =  z  to  obtain  the 
boundary  condition, 


d<(>  .  . 

V  =  ^  ^continuous 
z  az 


or  using  Equation  (86)  gives 


du 

r 

dz 


=  continuous. 


(104) 


(105) 


Continuity  of  the  Pressure 

Assuming  that  there  is  no  source  in  the  infinitesimal  volume  of  this 
cylinder,  the  volume  integration  of  Equation  (12)  gives 


^  p  n  dA  =  p  ^  J  V  dA  Ax 


(106) 


and  letting  Ax  — >  0  yields  the  pressure  boundary  condition 


Plls=  Pz's 


(107) 


which  means  that  the  acoustic  pressure  at  the  interface  is  continuous. 
Substituting  Equation  (86)  into  Equation  (28)  and  this  one  into  Equation  (107) 
gives 


p  u  =  continuous. 

n 


(108) 
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BOUNDARY  BETWEEN  LIQUID  AND  VACUUM 

As  a  very  good  approximation,  we  may  consider  this  boundary  as  pressure- 
release  for  acoustic  waves  in  the  liquid  layer.  Therefore,  the  only  boundary 
condition  is  that  the  pressure  vanishes  at  this  boundary.  Substituting 
Equation  (28)  for  the  pressure,  and  Equation  (86)  for  the  scalar  potential 
gives 


u  I  =  0.  (109) 

n  '  2 

0 

The  same  boundary  conditions  are  obtained  using  the  particle  velocity  and 
the  stress  tensor  of  an  elastic  layer  by  making  the  shear  speed  vanish.  The 
particle  velocity  is  given  from  Equations  (55)  and  (60).  From  Equations  (42) 
the  strain  tensor  in  an  axially  symmetric  environment  is  given  by 


a? 


r 


5r 


0 


S  = 


0 


0 


1 

2 


3?  5? 

r  2 

dz  *  dr 


0 


(110) 


and  using  Equation  (45)  for  the  stress  tensor  in  terms  of  the  components  of 
the  particle  velocity  gives 
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iw  ?  = 


5v 

(A+Z/i)  ^  ^ 


f  y  3v 
I  r  z 

[  F  ^ 


/•  3v  5v 


u  V  uv  \ 

ar'ar] 


/  3v  V 

dV  X 

0 

x\  ^  +  — 

+  ]  0 

(  3r  r 

5z  J 

5v  x 

dv  . 

-  V  5v  X 

Z 

0 

(A.Zp)  ^  .  a[ 

r  r  1 

.  ~  ^  ^  J 

(111) 

Hermit ian 

matrix 

as 

expected  for 

being  an  observable 

quantity. 

Now 

we  can  calculate  the  boundary  condition  for  some  other  cases  involving  the 
solid  layers. 


BOUNDARY  BETUEEN  LIQUID  AND  SOLID 

In  the  case  of  liquid-solid  boundaries  we  have  three  boundary  conditions 
to  satisfy. 


Continuity  of  the  Normal  Particle  Velocity 


Using  the  particle  velocity  vector  in  Equations  (60)  and  (55)  takes  us  to 
the  equation 


aF 


(112) 


which  by  substituting  Equations  (89a),  (91),  and  (94)  gives 


dz 

where  z^  stands  for  the  depth  of  the  boundary  to  be  matched,  and  these 
functions  must  be  evaluated  at  this  position.  The  subscript  "1"  stands  for  the 
depth  function  in  the  liquid  layer  and  "s"  labels  the  depth  function  in  the 


dz 


,  2 

+  k  V 
z  n  n 


(113) 
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solid  one. 


Continuity  of  the  Component  of  Stress 

Extracting  the  P  component  of  the  stress  matrix  in  both  media  gives 


i,p  (j)  <p  \  = 

0 


Lu 


dv 

(X  +  2»i)  ^  \ 


(114) 


where  by  the  same  equations  as  used  before  we  obtain 


P  u 

1  nl 


=  p  u 

z  s  ns 
0 


-  2p  (k  /k  )‘ 
z  s  n  s 

O 


u  + 

ns 


(115) 


The  P  Component  of  Stress  Vanishes 
rz 

From  Equations  (111),  and  the  same  equations  used  before,  this  boundary 
condition  becomes 


du 


dz 


(2k^  -  k^) 

n  s 


=  0. 


(116) 


BOUNDARY  BETWEEN  VACUUM  AND  SOLID 

For  a  free  elastic  layer  we  have  two  boundary  conditions  to  satisfy. 
These  are. 
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The  P  CorriDonent  of  Stress  Vanishes 
22 


u 

ns 


2 

0 


-  2(k  /k  f 

n  6 


( 


u  + 

ns 


(117) 


The  P  Comoonent  of  Stress  Vanishes 
rz 


2 


du 

ns 

dz 


+  (2k^ 

2  n 


0 


k^)  V 
s  n 


=  0. 


(118) 


BOUNDARY  BETWEEN  SOLID  LAYERS 

In  this  case  we  have  two  boundary  conditions  from  the  particle  velocity 
and  two  from  the  stress.  These  are  the  roost  general  of  all  boundary 
conditions  given  above  and  are  stated  as  follows: 


Continuity  of  the  Normal  Particle  Velocity 

du  ^ 

— +  k  V  =  continuous  (119) 

dz  n  n 


Continuity  of  the  Tangential  Particle  Velocity 

dv 

u  +  =  continuous  (120) 

n  QZ 
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where  it  is  emphasized  that  the  axially  symmetric  cylindrical  coordinate 
system  has  been  used  and  the  equations  are  implicitly  evaluated  at  the  depth 
of  the  boundary  interface. 
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CHAPTER  6 

THE  PROPAGATION  AND  MATCHING  ALGORITHM 


The  eigenvalues  that  characterize  the  stratified  ocean  are  obtained  by 
satisfying  all  the  boundary  conditions  simultaneously.  The  algorithm  developed 
here  consist  in  multiplying  all  the  calculated  propagation  and  matching 
matrices  to  obtain  the  value  of  the  compressional  and  shear  eigenfunctions  at 
some  interface  where  the  chosen  characteristic  equation  is  to  be  evaluated.  A 
combination  of  the  up-layer  and  down-layer  evaluation  of  the  eigenfunctions 
and  their  derivatives  will  be  used.  At  the  interface  between  two  solid 
layers,  the  four  boundary  conditions  in  Equations  (119)  through  (122)  are 
written  in  matrix  form  as 


J-i 


(Z  )  1 

(z  ) 

nJ-1 

3 

n3 

3 

u' 

(z  ) 

u' 

(z  ) 

nj-1 

3 

=  B 

3 

n) 

3 

V 

(z  ) 

V 

(z  ) 

nJ-1 

3 

n3 

3 

v' 

(z  ) 

v' 

(z  ) 

nJ-1 

J  J 

1  "J 

J 

(123) 


where 


B  = 
J 


1 

0 


P  T 
J  J 


0 

1 

2P. 


-p  T  /K 
3  3  3 


1 

0 


-P  Q 
3  3 


(124) 


and  it  is  designated  that  Q  s  2  k  /  k  and  T  s  1  -  Q  . 

3  n  3  3  3 


matching  of  the  eigenfunctions  the  matrix  B 


3-1 


For  up-layer 
in  Equation  (123)  is  taken  to 


the  other  side  and  multiplied  by  B^  to  obtain  the  matching  matrix 


38 


NSWC  TR  89-170 


where 


M  5  B  ^  B 
}  J-i  J 


M 

1 1 

0 

0 

M 

14 

0 

M 

22 

M 

23 

0 

0 

M 

32 

M 

33 

0 

M 

41 

0 

0 

M 

4  4-' 

M=M=Q  +TR, 

11  33  J-1  J  J 


M  =  Q  -  Q  R  , 

14  j-i  J  j’ 


M=M=T  +Q  S. 

22  44  J-1  J-1  J 


(125) 

(126a) 

(126b) 

(126c) 


and 


M  =  T  -  T  R  , 

41  J-1  J  J 


M  =  M  k  . 

23  41  n 


M  =  2  (1  -  S  )/K 

32  J  J-1 


(126d) 

(126e) 

(126f) 


are  the  elements,  and  the  ratio  of  densities  and  shear  moduli  are  defined  as 
R  =  p  /p  and  S  =  u  /u 

J  ^j-i  J  ^j  "^j-i 

At  the  liquid-liquid  interface  the  boundary  conditions  are  reduced  to 
two,  and  the  up-layer  matching  matrix  becomes 


W  = 
J 


0  ' 
1 , 


(127) 


at  the  interface. 

These  matching  matrices  simply  provide  the  eigenfunctions  and  their 
derivatives  at  the  bottom  of  the  J-1  layer  when  the  values  at  the  top  of  the 
layer  are  known.  Another  matrix  is  needed  to  propagate  the  solution  from 
the  bottom  to  the  top  of  the  layer. 

When  the  values  at  the  bottom  of  the  layer  has  been  determined,  the 
solution  of  the  wave  equations  can  be  used  to  propagate  the  solution  to  the 
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top  of  this  layer.  To  create  this  propagation  matrix,  Equations  (95)  and  (100) 
are  rewritten  in  the  matrix  form 


fu  (z) 
nj 

U'  (Z) 
nJ 

V  (z) 
nJ 

V'  (z) 
nJ 


^  Ai[C^j(z)]  Bi[<^^(z)]  0 

-S^^^Ai' [c  (z)]  (z)]  0 

J  nJ  J  nJ 

0  0  expdy  ) 

nJ 

0  0  iy  expdy  ) 

n  j  n  J 


0 

fa 

nJ 

0 

b 

nJ 

exp(-fr  ) 
nJ 

C 

nJ 

iy  exp(-i3r  ) 
n  J  nJ  J 

d 

nJ 

(128) 


where  the  banded  4x4  matrix  means  that  the  shear  and  compressional 
eigenfunctions  are  propagated  independently  from  each  other  even  though  they 
depend  on  each  other  in  the  matching  process.  Therefore,  the  propagation 
matrix  can  be  divided  into  a  2x2  compressional  propagation  matrix  and  a  2x2 
shear  propagation  matrix.  Inverting  and  evaluating  the  compressional 
propagation  matrix  at  the  bottom  of  the  j*'*’  layer  gives 


n  Bi'  (< 
-n  Ai'  [<; 


nJ 

nJ 


(z 

J 

(z 


J+l 

'j+1 


)] 

)] 


-71  S-'^^Ail^  (z  )] 
J  nJ  J+1 


U' 

nJ 


)' 

) 


(129) 


where  the  Wronskian  relationship 


54 


Bi'((:;)  Ai(<)  -  Bi(<)  Ai'(C)  =  I/ti 


(130) 


has  been  used  and  the  primes  denote  the  derivative  with  respect  to  the 
argument.  For  the  shear  propagation  matrix  we  obtain 


fc 

'1/2  exp(-iy  z  )  -i  exp(-ty  z  )/{2y  )' 

'v  (z 

nJ 

nJJ+l  nJJ+1  nJ 

nJ 

J*1 

d 

njJ 

1/2  exp(iy  z  )  i  expily  z  )/(2y  ) 

[  nJ  J+l  ^  nj  J+1  nj  J 

v'  (z  ) 

[  nJ  j 

(131 ) 
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where  substituting  the  unknown  coefficients  back  into  Equation  (128)  and 
evaluating  it  at  the  top  of  the  layer  gives 


fu  (z  )  1 
nJ  J 

c 

11 

c 

12 

0 

0  ' 

u'  (z  ) 
nJ  J 

c 

21 

c 

22 

0 

0 

V  (z  ) 
nJ  J 

0 

0 

c 

33 

c 

34 

< 

N 

4^ 

0 

0 

c 

43 

c 

44 

/ 

(z 

nJ 

J  +  l 

u' 

(z 

nJ 

J+l 

V 

(z 

nJ 

J+l 

v' 

(z 

1  "J 

J+l 

(132) 


where  the  elements  of  the  compressional  eigenfunction  are  given  by 


C  = 
n 


[  Ai[^  (z  )] 

I-  nj  J 


Bi'  (C  (z.  )]  -  Ai' 

nj  J+ 1 


(2  )] 
nJ  j+1 


Bi[C  (z  )] 
nj  J 


(133a) 


C  =  n  r  Ai(<  (z  ))  Bi(C  (z  ))  -  Ai(<  (z  ))  Bi(C  (z  ))  ]/S 

12  L  ^nj  J  nj  j+1  nJ  J*1  ^nj  J  -•  J 


1/3 


(133b) 


C  =  ir  r  Ai'(C  (Z  ))  Bi'«  (z  ))  -  Ai' (C  (z  ))  Bi'(<:  (z  ))  ]S 

21  L  ^nj  J+1  ^nj  J  ^nj  J  ^nj  J+1  J  J 


1/3 


(133c) 


and 


C  =  n  r  Ai[<  (z  )]  Bi' [<  (z  )] 

22  L  ^nj  J+l  nJ  J 


(133d) 


and  the  elements  for  the  shear  eigenfunction  are  given  by 


C  =  cosij  D  ) 

44  nj  J 


(134a) 


and 


C  =  -y  '  sin(y  D  ) 

34  nJ  nJ  J 

C  =  y  sin(y  D  ) 

43  nJ  nJ  J 


(134b) 

(134c) 


where  D  =  z  -  z  is  the  thickness  of  the  layer  and  y  is  given  after 

J  j  +  i  J  ^  nj  ® 

Equation  (100).  For  trapped  compressional  modes,  only  the  exponentially 
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decaying  solution  of  the  Airy  functions  must  be  used  since  there  is  no 
acoustic  source  in  the  bottom.  If  Ai(^)  is  evaluated  at  the  exponentially 
decaying  part  of  the  function,  then  Bi(C)  cannot  be  evaluated  simultaneously 
because  it  is  exponentially  increasing  producing  a  floating  point  overflow  in 
the  computer.  Therefore,  an  Airy  function  subroutine  that  can  calculate  one  of 
the  solutions,  instead  of  both  of  them  simultaneously,  must  be  used  Another 
serious  complication  that  occurs  with  the  evaluation  of  Equations  (133)  is  the 
numerically  unstable  result  when  the  subtraction  of  large  but  very  close 
numbers  is  performed.  In  the  case  where  |C|  >>  10  and  Arg  <  =  n  the  Airy 
functions  have  the  asymptotic  forms. 


Bi(<)  L  Ai(C)  (135a) 

and 

Bi' (<)  — )  i  Ai' (<)  (135b) 

which  causes  precision  problems  when  evaluating  the  Wronskian  and  the  elements 
of  the  propagation  matrix.  These  same  propagation  matrix  coefficients  have 
been  encountered  by  Gordon^^’^*  in  the  area  of  quantum  mechanics  and  he  solved 
the  problem  by  direct  substitution  of  the  series  solutions  for  the  Airy 
functions  and  its  derivatives  in  Equations  (133)  to  factor  out  the  exponential 
or  sinusoidal  functions.  However,  these  where  performed  under  the  assumption 
of  purely  real  arguments  of  the  Airy  functions  and  the  series  solutions  used 
do  not  apply  to  complex  arguments.  A  similar  substitution  and  cancellation 
technique  is  obtained  by  evaluating  the  numerically  stable  independent  Airy 

57 

results  given  by  Schulten,  Anderson,  and  Gordon  to  avoid  floating  point 
overflow  and  other  precision  problems. 

In  the  water  layers  the  propagation  relationship  reduces  to 


f  u  (z  )1 

fc 

c  1 

'  u  (z  )' 

nj  J 

11 

12 

nJ  J+1 

u'  (z  ) 

c 

c 

U'  (z  ) 

i  nJ  J  J 

21 

22j 

4- 

c 

where  the  unknowns  are  given  in  Equations  (133). 

In  the  models  to  be  described  next,  it  is  further  assumed  that  the 
acoustic  properties  of  each  elastic  layer  are  constant.  This  is  done  for 
simplicity  and  with  the  knowledge  that  it  is  not  easy  to  determine  all  the 
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fine-structure  properties  of  the  bottom  sediments  experimentally,  therefore 
the  bottom  is  usually  described  as  layers  of  constant  properties^^’^^ 


CASE  A.  RIGID  FALSE  BOTTOM  MODEL 

The  rigid  bottom  interface  is  at  j=F  in  Figure  4,  the  water-bottom 
interface  is  at  j=J,  and  the  soft  surface  of  the  ocean  is  at  J=l.  The  r5gid 
"false"  bottom  interface  is  artificially  incorporated  only  to  convert  a 
continuous  wave  number  spectrum  into  a  discrete  form.  Therefore,  the  trapped 
modes  will  be  treated  as  if  the  rigid  bottom  does  not  exist.  Each  water  layer 
has  two  unknowns  to  be  determined,  and  each  solid  layer  has  four.  The  rigid 
bottom  has  null  tangential  and  normal  particle  velocity.  Therefore,  the 
boundary  conditions  for  the  radiating  modes  at  the  "false"  bottom  interface 
are 


and 


u  (z  ) 

nF-l  F 


v'  (z  )  =  0 
nF-l  F 


(137a) 


u'  (z  ) 

nF-l  F 


+  E  V  (z  )  =  0 
n  nF-l  F 


(137b) 


2 

which  have  four  unknowns  for  any  trial  value  H  of  the  eigenvalue  squared  k  . 

n  n 

These  equations  are  conveniently  written  in  the  matrix  form 


fu  (z  ) 

'  0  -1  ' 

U. 

1 

U. 

c 

u'  (z  ) 

-E  0 

nF-l  F 

— 

n 

V  (z  ) 

1  0 

nF-l  F 

V'  (z  ) 

0  1 

u. 

1 

u. 

c 

.  ! 

V  (z  ) 
nF-l  F 

V'  (z  ) 
nF-l  F 


(138) 


where  the  special  4x2  matching  matrix  is  represented  by  W^. 

This  set  of  equations  is  substituted  into  the  up-layer  propagation  matrix 
for  the  F-1  layer  to  obtain 
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[U  (z  ) 
nF-l  F-1 


/  f  s  V  (Z  ) 

U  (Z  )  ,r>  lu  nF-l  F 

nF-l  F-1  =  C  W  ,  r  ^ 

,  V  F-1  F  V'  (Z  ) 

V  (z  )  nF-l  F 

nF-l  F-1 

V'  (Z  ) 

nF-l  F-1 


(139) 


where  both  eigenfunctions  and  their  derivatives  are  still  to  be  determined. 
This  new  recurrence  relationship  is  substituted  into  the  up-layer  matching 
matrix  for  the  F-1  solid-solid  boundary  to  obtain  the  new  relationship 


fu  (z  ) 

nF-2  F-1 

U'  (Z  ) 

nF-2  F-1 

V  (z  ) 

nF-2  F-1 

V'  (Z  ) 

nF-2  F-1 


V  (z  ) 

=  W  CM 

F-1  F-1  F  V'  (Z  ) 
nF-l  F 


(140) 


and  by  the  same  token,  the  matching  and  propagating  matrices  for  the  solid 
layers  are  multiplied  to  each  other  until  the  water-bottom  interface  is 
reached  with  the  relationship 


=  c  w  c 

J  J+1  J+1 


W  C  M 

F-1  F-1  F 


V  (Z  )) 
nF-l  F 

V'  (Z  )  ■ 

nF-l  F 


fu  (Z  ) 
nJ  J 

u',(z  ) 

n  J  J 

V  (z  ) 
nJ  J 

V'  (z  ) 
n  J  J 


The  multiplication  of  the  propagating  and  matching  matrices  in  the  solid 
layers  is  now  given  by  the  4x2  matrix 


(141) 


EsCW  C  ...  Vi  C  W. 

J  J+l  J+l  F-1  F-1  F 


(142) 


Next  is  a  use  of  the  liquid-solid  boundary  conditions  to  evaluate  the 
compressional  eigenfunction  at  the  water  column.  The  three  boundary  conditions 


u'  (z  )  =  u'  (z  )  +  E  V  (z  ) 

nJ-l  J  nJ  J  n  nJ  J 


(143a) 


p  u  (z)=p  T  u  (z)-p  Q  v'  (z)  (143b) 

'^J-l  nJ-l  J  J  nJ  J  J  nJ  J 
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and 


2  u'  (z  )  +  (2  E  -  ic  )  V  (z  )  =  0 
nJ  J  n  J  nJ  J 


(143c) 


where  the  first  two  equations  are  needed  to  evaluate  the  compressional 
eigenfunction  in  the  liquid  layer,  and  the  third  condition  will  be  used  as  the 
characteristic  equation  of  the  environment. 

The  first  two  equations  are  rewritten  in  matrix  form  as 


u  (z  )' 

nJ-1  J 

r  T  /R 

J  J 

0 

0  -Q  /R  1 

J  J 

U  '  (z  ) 
nJ-1  J 

0 

1 

E  0 

u 

(z  ) 

n  J 

J 

u' 

(z  ) 

nJ 

J 

V 

(z  ) 

n  J 

J 

v' 

(z  ) 

n  J 

J 

(144) 


where  the  2x4  special  matching  matrix  will  be  known  as  Substituting 

Equations  (137)  into  this  set  of  equations  gives 


nJ-l  J 

u  '  (z  ) 

nJ-1  J 


fv  (z  )] 

=  W  E 

nF-1  F 

V'  (z  ) 

J 

u. 

1 

u. 

c 

(145) 


Now  comes  the  propagation  and  matching  in  the  liquid  layers.  The  method 
is  the  same  as  for  the  solid  layers,  except  that  now  the  matrices  to  multiply 
are  2x2.  After  propagating  and  matching  up  to  the  surface  it  is  found  that 


r  u  (z  )i 

V  (z  )' 

nl  1 

=  C  M  C  M  C  .  . 

,  .  M  C  ME 

nF-l  F 

v'  (z  ) 

1  ">^-1  ^ 

u'  (z  ) 

nl  1 

V  / 

1  2  2  3  3 

1 

1 

(146) 


where  the  final  product  of  all  the  propagation  and  matching  matrices  is 
defined  as  the  2x2  matrix 


FSCIMC1MC...1M  C  ME.  (147) 

1  2  2  3  3  J-1  J-1  J 

The  pressure-release  boundary  condition  of  the  surface  (J=l)  gives 

u  (z  =  0)  =  0  and  the  derivative  of  the  eigenfunction  at  the  surface  will  be 
n  1  1 

arbitrarily  set  to  unity  since  the  normalization  condition  will  take  care  of 
its  proper  evaluation.  This  gives 


45 


NSWC  TR  89-170 


'  0 
.  1 


=  F 


V 

nF-l 


v' 

nF-l 


(148) 


where  inverting  F  gives 


and 


V  (z  ) 

nF-l  F 


(F  -  F  (F  /F  )) 
22  12  21  11 


-1 


V'  (Z  ) 
nF-l  F 


(F  -  F  (F  /F  )) 
21  11  22  12 


-1 


(149a) 

(149b) 


A  final  substitution  of  the  shear  eigenfunctions  at  the  rigid  bottom 
interface  into  Equation  (141)  gives  the  two  necessary  values  to  be  substituted 
into  the  chosen  characteristic  equation.  Equation  (143c).  The  trial  value  E 

n 

is  the  square  of  an  eigenvalue  k  when  the  complex  characteristic  equation  (or 

n 

determinant)  is  null,  i.e.. 


Wik  )  £  2  u'  (z  )  +  (2  -  K^)  V  (z,)  =  0.  (150) 

n  nJ  J  n  J  nJ  J 

Note  that  the  objective  is  to  calculate  only  u' ,(z,)  and  v  ,(z,).  Therefore, 
there  is  no  need  to  calculate  the  eigenfunctions  at  the  other  interfaces  to 
find  the  eigenvalues  of  the  problem.  This  same  method  of  matrix  multiplication 
is  also  used  to  evaluate  the  eigenfunctions  at  all  the  interfaces  for  the 
normalization  calculation  described  in  the  next  chapter. 


CASE  B.  SEMI -INFINITE  BASEMENT  MODEL 

This  model  is  used  for  the  trapped  modes  and  the  radiating  modes 
undiscriminatingly.  The  trapped  modes  have  the  property  that  they  are 
exponentially  decreasing  with  depth  in  an  isovelocity  layer  and  the  radiating 
modes  represent  out-going  propagating  waves  that  oscillate  towards  infinity 
without  reflections  but  with  absorption  which  causes  its  damping.  Both, 
compressional  and  shear,  eigenfunctions  can  be  trapped  or  radiating. 

Therefore,  it  is  no  longer  appropriate  to  refer  to  a  mode  as  simply  trapped  or 
radiating.  In  fact,  there  are  four  mode  classifications.  There  is  the  most 
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common  mode  where  the  compressional  eigenfunction  is  trapped  and  the  shear 
eigenfunction  is  radiating.  This  type  of  mode  will  be  labeled  a  "T-R  mode," 
where  the  first  letter  always  refers  to  the  compressional  function.  The 
compressional  eigenfunction  of  the  T-R  mode  in  the  semi-infinite  layer  is 
described  by  the  exponential  function 

u  ^  (z)  =  a  exp(-7)  ^  z)  (151a) 

nF- 1  nF“ 1  nF-1 

and 

u'  (z)  =  -7)  ^  ,  u  ^  ,  (z)  (151b) 

nF-1  nF-1  nF- 1 

2  2  2 

where  t)  =  k  -  k  .  The  shear  eigenfunction  is  radiating  without 

nF-1  n  F-1 

reflections,  i . e. , 


and 


V 

nF-  1 


(Z) 


c  exp(i9f  z) 

nF-1  ^  ®nF-l 


V' 

nF-1 


(z) 


=  tr 


nF-1 


V 

nF-1 


(z)  . 


(151c) 

(151d) 


Now  the  eigenfunctions  at  the  top  of  the  basement  layer  are  given  by 


'  u  (z  )  ' 

nF-1  F-l 

'  1 

0 

u'  (z  ) 

nF-1  F-l 

= 

"V-l 

0 

V  (z  ) 

nF-1  F-l 

0 

1 

v'  (z  ) 

nF-1  F-l  J 

0 

iy 

nF-  l-* 

V 

u  (z  ) 

nF-1  F-l 

V  (z  ) 

nF-1  F-l 


(152) 


which  gives  the  values  at  the  top  of  the  basement  layer.  Lets  define  the 
special  4x2  propagating  matrix  in  Equation  (152)  as  ^  and  only  two  of  its 
elements  must  be  changed  if  the  mode  becomes  a  T-T,  R-T,  or  R-R  mode.  To 
obtain  the  values  at  the  bottom  of  the  F-2  layer,  we  multiply  by  the  matching 
matrix  at  this  interface  which  gives 


u  (z  )  1 

nF-2  F-l 
u'  (z  ) 

nF-2  F-l 

V  (z  ) 

nF-2  F-l 

V'  (Z  ) 

nF-2  F-l 


(M  C 

F-l  F-l 


1 

(z 

nF-1 

F-l 

(z 

nF-1 

F-l 

(153) 
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which  corresponds  to  Equation  (140)  in  the  previous  matrix  multiplication 
method.  By  the  same  matching  and  propagating  algorithm  the  values  at  the 
surface  of  the  ocean  wave  guide  are  reached  with  the  relationship 


'  0  ' 

u 

(z  )  1 

=  F 

nF-l 

F-1 

(154) 

1 

V  J 

i  V  (z  )  J 

^  nF-l  F-1  ^ 

where  the  2x2  matrix 


F=C[MC...IMC...IM  C  (155) 

12  2  J  J  F-1  F-1 

is  inverted  to  obtain  u  (z  )  and  v  (z  ).  After  calculating  their 
derivatives  using  Equations  (151b)  and  (151d),  the  solution  is  propagated  up 
to  the  top  of  the  layer  to  obtain  u'  (z  )  and  v  (z  )  which  are  n^^eded  to 

nJ  J  nJ  J 

calculate  the  characteristic  equation,  Equation  (150).  This  method  requires 
one  loop  of  matrix  multiplications  in  the  solid  to  obtain  the  matrix  E  in 
Equation  (142),  which  is  saved  for  later  use,  but  is  also  used  to  keep 
multiplying  matrices  in  the  liquid  layers. 

After  evaluating  the  eigenvalue,  the  eigenfunction  may  be  calculated  by 
down-layer  multiplication  of  propagation  and  matching  matrices.  For  this 
purpose,  the  matrices  are  slightly  different  from  the  up-layer  ones.  The 
elements  of  the  matching  matrix  in  Equation  (125)  for  down-layer  matching  of 
the  eigenfunctions  are  given  by 


J,  4, 

M  =M  =Q  +T/R  , 

n  33  j+i  J  j*i 


M  =  Q  -  Q  /  R  , 

14  J+1  J  J  +  1 


M  =M  =T  +Q  /S  , 

22  44  J+1  J+1  J+1 


M  =  T  -  T  /  R  , 

41  J+1  J  J+1 


(156a) 

(156b) 

(156c) 

(156d) 
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and 


4^  ^  p 

M  =  M  k  . 

23  41  n 


=  2  (1  -  1  /  S  )  /  , 

32  J*1  J+1 


(156e) 


(156f) 


The  arrows  have  been  include  to  emphasize  the  down-layer  propagation  and 
matching  matrices  and  to  distinguish  them  from  the  up-layer  ones.  The 
compressional  elements  of  the  down-layer  propagating  matrix  in  Equation  (132) 
are  given  by 


C  =  n  [  AifC  ,(2,^  Bi'[<  (2  )]  -  AiMC  )]  Bi[<  (z  )]  ] 

11  nj  J  +  1  nj  J  nj  J  nj  J  +  1  J 


(157a) 


C,  =  IT  [  Ai«  (z  ))  Bi(C  (z  ))  -  Ai(C  (z  ))  Bi(C  (z  )) 

12  L  J+1  '’nj  j  '’nj  j  ^nj  J-l  ■>  J  ’ 


(157b) 


K  ■ 

(157c) 


and 


C  =  n[Ai[<(z)]Bi'[^(z  )]-Ai'[C(z  )]Bi[C  (z)]] 

22  ^  nJ  J  nJ  J+1  nJ  J+1  ^nj  J  J 


The  shear  elements  of  the  down-layer  propagating  matrix  are  given  by 


(157d) 


4^  •4- 


and 


=  cos(y  D  ) 

44  nJ  J 

(158a) 

*  sin(y  D  ) 

(158b) 

nJ  nJ  J 

7  sin(9r  D  ). 

(158c) 

However,  to  obtain  the  eigenfunctions,  the  zeros  of  the  complex 
characteristic  equation.  Equation  (150),  must  be  found.  The  most  optimum 

4  6 

method  for  finding  the  complex  zeros  may  never  be  known  since  many 
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challenging  environmental  conditions  can  be  encountered.  Examples  are 
secondary  sound  channels  and  surface  ducts  which  cause  degenerate  eigenvalues 
with  irregular  spacing  and  some  of  these  eigenvalues  may  be  too  close  to  each 
other  for  the  limited  precision  of  a  computer.  Ellis^°  has  tackled  a  similar 
problem  with  a  two-ended  shooting  technique,  but  the  model  incorporates  the 
effects  of  shear  wave  and  the  attenuation  coefficients  as  an  approximation  for 
small  values,  he  does  not  take  into  account  the  radiating  modes,  and  does  not 
search  for  the  eigenvalues  in  the  complex  plane.  However,  a  simple  algorithm 
may  be  obtained  after  analyzing  the  behavior  of  Equation  (150)  for  a 
simplistic  ocean  environment  and  it  may  be  eventually  improved  as  we  encounter 
more  difficult  situations. 

The  complex  determinant  in  Equation  (150)  was  calculated  as  a  function  of 

the  complex  trail  wave  number  for  the  very  simple  case  of  a  water  column  over 

a  semi-infinite  bottom  layer.  To  start  with  a  very  simple  case,  the  bottom  is 

a  fluid  layer  with  the  acoustic  properties  in  the  second  row  of  Table  1.  These 

1  7 

fluid  properties  are  obtained  from  Jensen,  the  properties  of  clay-silt,  sand, 
and  basalt  are  taken  from  the  paper  by  Werby  and  Tangof'  and  the  parameters 
for  chalk  are  given  by  Ellis! ^  The  first  case  is  the  fluid  bottom  with  a  water 
column  200.0  meters  deep  and  a  constant  sound  speed  of  1500.0  m/s.  The  bottom 
is  fluid-like  because  of  the  very  small  shear  speed  (1.0  m/s)  and  the 
relatively  large  shear  attenuation  coefficient  to  rapidly  damp-out  the 
remaining  shear  contribution.  The  frequency  of  the  sound  is  25.0  Hz  and  all 
other  bottom  properties  are  given  in  Table  1. 

The  contour  plot  of  the  complex  determinant  is  displayed  in  Figure  6.  The 
dash  curves  represent  the  contour  where  the  real  part  of  the  determinant  is 
zero,  and  the  solid  curves  represent  the  contour  where  the  imaginary  part  of 
the  determinant  is  zero.  Since  both  components  of  Equation  (150)  must  be  zero 
simultaneously,  then  the  complex  eigenvalues  are  given  by  the  points  where 
both  curves  intersect.  Note  that  these  curves  are  perpendicular  at  the  point 
of  intersection,  therefore  Equation  (150)  is  an  analytic  equation  that 
satisfies  Cauchy-Riemann  relations.  This  simplifies  the  Newton-Raphson  method 
for  converging  into  the  complex  zeros.  The  real  part  of  the  determinant  along 
the  real  axis  is  plotted  in  Figure  7  where  the  number  of  eigenvalues  is  equal 
to  the  number  of  zeros  of  this  curve.  The  maximum  trail  wave  number  is  given 
by  the  minimum  speed  in  the  water  column  and  a  general  behavior  of  this  curve 
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is  that  it  tends  to  infinity  as  the  wave  number  approaches  zero.  A  search  for 
the  change  in  the  sign  of  this  curve  will  be  used  to  obtain  an  initial 

4  1 

estimate  of  the  eigenvalue.  Nagl  made  a  detailed  study  of  the  behavior  of 
the  purely  real  eigenvalues  along  the  real  axis  and  derived  the  useful 
approximation 


Ak(n)  s  (n  R  c  )/(2  f  -z)  (159) 

min  j 

which  gives  the  approximate  spacing  of  the  eigenvalues  as  a  function  of  the 
mode  number  n,  the  frequency  f,  the  depth  of  the  ocean  floor,  and  the  minimum 
sound  speed  in  the  water  column.  This  equation  is  used  to  establish  the  step 
size  in  the  search  for  the  first  couple  of  eigenvalues.  The  spacing  between 
the  previous  two  eigenvalues  is  later  used  to  establish  the  step  size  for  the 
next  eigenvalues. 

After  these  zeros  are  found,  we  vary  the  imaginary  part  of  the  trail  wave 
number  in  order  to  follow  the  dash  curves  in  Figure  6  until  a  change  in  sign 
of  the  imaginary  part  of  the  determinant  is  detected.  Newton-Raphson*^  method 
was  tried  to  pin-down  the  complex  eigenvalues  in  double  precision  accuracy  but 
there  where  times  when  it  converged  to  unrealistic  values.  This  method  was 
also  used  by  Otsubo  for  relatively  simple  ocean  environments  without  the 
effects  of  shear  waves  but  we  have  discarded  it  as  "not  uniformly  convergent." 
A  more  uniformly  convergent  method  is  described  by  Morris^®  where,  instead  of 
searching  for  the  complex  zeros  of  Equation  (150),  we  search  for  the  local 
minima  of  the  magnitude  of  the  determinant.  These  minima  are  located  exactly 
where  the  complex  eigenvalues  are,  the  method  uniformly  converges  to  the  local 
minima,  but  the  convergence  is  slightly  slower  than  the  Newton-Raphson 
method.  Since  the  uniform  convergence  is  more  important,  this  is  the  method 
used  in  this  model. 

A  plot  of  the  eigenvalues  in  the  complex  k-plane  is  given  in  Figure  8 
where  the  three  trapped  modes  have  eigenvalues  with  negligible  imaginary 
parts,  and  the  four  radiating  modes  have  eigenvalues  with  real  parts  smaller 
than  0.09  1/m.  Note  that  the  radiating  eigenvalues  have  such  a  large  imaginary 
part  that  they  get  rapidly  attenuated  as  they  propagate  in  range  therefore 
contributing  only  to  the  intensity  near  the  source. 

If  a  rigid  false  boundary  is  placed  at  a  depth  of  1200.0  meters?^  the 
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radiating  spectrum  becomes  the  one  displayed  in  Figure  9,  which  agrees  with 
Miller’s  perturbative  approximation.  The  real  part  of  the  minima  corresponds 
to  the  four  radiating  eigenvalues  in  Figure  8,  which  means  that  the  rigid 
false  bottom  method  is  trying  to  give  preference  to  the  true  radiating 
eigenvalues  over  the  "pseudo-eigenvalues"  by  giving  the  latter  a  larger 
imaginary  part. 

The  trapped  and  radiating  eigenvalues  are  not  unique  in  the  complex  wave 

number  spectrum.  Surface  waves  may  also  be  incorporated  in  this  ocean  model. 

The  surface  wave  produced  at  a  solid-vacuum  interface  is  called  a  Rayleigh 

wave.  In  the  case  of  a  liquid-solid  interface  the  surface  wave  is  called  a 

Generalized  Rayleigh  wave,  while  the  solid-solid  surface  wave  is  called  a 

Stoneley  wave.  These  waves  satisfy  the  wave  equation  and  the  boundary 

conditions  but  they  decay  exponentially  in  both  directions  very  rapidly. 

Therefore  surface  waves  are  not  expected  to  be  of  extremely  high  importance 

when  the  receiver  is  greater  than  a  wavelength  away  from  interface,  even 

though  they  may  be  added  in  these  transmission  loss  calculations  without 
,,  ..  63,64 

complications. 

The  spherical  waves  diverging  from  the  source  may  be  expressed  as  the 
infinite  summation  of  plane  waves.  Each  plane  wave  hits  the  liquid-solid 
boundary  and  each  solid-solid  interface  of  the  wave  guide  at  a  different  angle 
of  incidence.  Some  of  them  will  have  an  angle  of  incidence  greater  or  equal  to 
the  compressional  critical  angle  and/or  the  shear  critical  angle  of  an  elastic 
layer.  Under  this  condition,  total  internal  reflection  occurs  and  the  wave 
propagates  parallel  to  the  interface.  This  wave  is  also  called  a  surface  wave 
since  its  amplitude  decays  exponentially  with  the  normal  distance  from  the 
interface  along  which  the  propagation  occurs. 

Lateral  waves  propagate  at  the  solid  side  of  the  liquid-solid  interface 
and  they  are  automatically  incorporated  in  the  model  when  the  trapped  and 
radiating  modes  are  used.  When  a  spherical  wave  from  a  nearby  source  hits  the 
liquid-solid  interface,  two  lateral  waves  are  created  which  will  propagate  in 
the  elastic  layer  and  parallel  to  the  interface.  One  is  caused  by  the  shear 
critical  angle  of  incidence  and  the  other  is  caused  by  the  compressional  onef^ 
As  they  propagate,  they  reradiate  back  into  the  liquid  layer  causing  the 
Schmidt  head  wave. 
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CHAPTER  7 

THE  NORMALIZATION  COEFFICIENT 


The  depth  functions  are  complete  and  orthogonal.  However,  they  are  not 
necessarily  normalized  since  the  wave  equation  is  satisfied  regardless  of  the 
constant  in  front  of  the  solution.  Remember  that  the  derivative  of  the 
eigenfunction  at  the  surface  has  been  set  to  unity  under  the  condition  that 
the  normalization  constant  would  take  care  of  this  last  unknown.  For  the 
transmission  loss  calculation  in  the  next  chapter,  the  summation  of  normalized 
eigenfunctions  is  needed  to  obtain  the  proper  contribution  of  each  mode. 

The  normalization  constant  for  each  mode  is  given  by 


N 

n 


F-1 

I 


j=i 


N 

nj 


(160a) 


where  N  is  the  contribution  of  the  j''*'  layer  to  the  n*'*'  mode,  i.e., 
nj 


N 

nj 


-Z 

=p  u  (z)u  (z)dz. 

J  nj  nj 

Z'' 

J 


(160b) 


The  normalization  coefficient  is  in  general  a  complex  number,  and  the 
eigenfunctions  at  all  the  interfaces  are  known  after  the  eigenvalues  are 
found. 

In  this  ocean  model,  the  water  layers  have  variable  sound  speed  and 
negligible  attenuation  coefficient  compared  to  the  absorption  of  the  bottom 
layers.  However,  these  elastic  layers  of  the  bottom  have  constant  acoustic 
properties.  Therefore,  the  normalization  calculation  for  both  types  of  layers 
is  different. 

The  water  layers  are  defined  to  be  the  ones  where  1  s  j  <  J.  Gordon’s 
formulas^^  provide  the  analytical  integration  of  a  linear  combination  of  Airy 


53 


NSWC  TR  89-170 


functions.  The  formulas  to  use  are 


and 


where 


and 


represent  any  linear  combination  of  Airy  functions. 

Equation  (161a)  relates  to  Equation  (160b)  where  A  =  B.  In  this  case, 


where 

&  H  )  /  S  -  z  .  (164) 

nj  J  n  J  J 

Note  that,  even  if  the  attenuation  coefficient  of  the  water  column  is 
neglected,  the  argument  of  the  Airy  functi'-n  is  complex  because  of  the  complex 
eigenvalue.  Therefore,  the  absorption  of  'he  bottom  layers  are  affecting  the 
propagation  of  sound  in  the  water  column 

The  calculation  of  the  normalization  contribution  in  the  elastic  layers 
is  simplified  by  the  assumption  of  sediments  with  constant  acoustic 
properties,  but  a  distinction  must  be  made  for  the  case  where  the  function 
decays  exponentially  with  depth  (trapped  mode)  or  if  it  oscillates  (radiating 
mode ) . 

In  the  case  of  a  trapped  mode,  the  compressional  eigenfunction  in  the 
layer  is  given  by 

u  (z)  =  u  (z  )  expf-T)  (z  -  z  )1  (165) 

nJ  nJ  J  nJ  J  J 


I 


A[a(z+b)]  B[a(z+b)]  dz  =  (z+b)AB  -  A’B’/a 


(161a) 


r  A[a(z+b  )]  B[a(z+b  )]  dz  = 

J  1  2  2 


A’B  -  AB’ 


(161b) 


a  (b  -  b  ) 
1  2 


A[a(z+b^ ) Isa^Ai [a(z+b^ ) ]+b^Bi [a(2+b^ ) ] 


(162a) 


B[a(z+b  )]=a  Ai[a(2+b  )]+b  Bila(2+b  )] 
2  2  2  2  2 


(162b) 
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2  2  2 

for  z  £  z  £  z  and  where  t)  =  k  -  k  .  Substituted  into  Equation  (160b) 

J  J  +  l  nj  n  J  V  ^ 

gives  the  solution 


N  =p  u^(z)il  -  exp  [-2  D  n  1 1  /  [  2  t?  ] 
nJ  nJ  J  \  J  njJJ  nJ 


(166) 


where  D  =  z  -  z  is  the  thickness  of  the  elastic  layer. 

J  J+i  J 

The  compressional  eigenfunction  to  use  for  the  oscillatory  depth  function 
in  the  layer  is 


u  (z)  =  u  (z  )  cos[7)  (z-z  )]  +  u'  (z  )  7)*^  sinLT}  (z-z  )]  (167) 

nJ  nJ  J  nJ  J  nJ  J  nJ  nJ  J 


and  the  Equation  (160b)  becomes 


(  2  ,2  -zUl  r.  .  ^  "nj'’ 

N  =  p  u  +  u  T)  -  D  +  - T - - - —\  +  p  u  u  - - - — 

nJ  nJ  nJ  njJ[2  J  4  7)  J  *^1  nj  nJ  7) 

^  \  ^  4 


'sin(D  7)  )]' 

J  nJ 


(168) 


where  the  argument  of  the  eigenfunction  has  been  omitted  for  simplicity  and 
where  J  £  J  <  F. 

Finally,  the  normalized  eigenfunctions  and  their  derivatives  are  given  by 


u  (z)  =  u  (z)  / 

nj  nj  n 


(169a) 


u'  (z)  =  u'  (z)  / 

n  J  n  J  n 


(169b) 


V  (z)  =  V  (z)  /  N 

n  J  nJ  n 


(169c) 


v'  (z)  H  v'  (z)  /  N 

n  J  nJ  n 


(169d) 


After  all  the  eigenvalues  and  normalized  eigenfunctions  are  found,  it  is 
next  to  calculate  the  transmission  loss  of  sound  in  this  ocean  wave  guide. 
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CHAPTER  8 

THE  TRANSMISSION  LOSS 


Any  massive  object  that  vibrates  radiates  acoustic  energy.  Power  is  the 
time  rate  at  which  energy  is  radiated  and  intensity  is  defined  as  the  rate  of 
energy  flow  through  a  unit  area.  The  intensity  is  a  vector  quantity  which  also 
gives  the  direction  of  the  energy  flow.  However,  the  fluctuations  of  the 
intensity,  the  pressure,  and  the  power  are  of  order  of  magnitudes  and  this 
presents  problems  when  plotting  them.  Therefore,  in  acoustics,  the  intensity 
is  converted  into  decibels  in  order  to  reduce  the  high  fluctuations.  The  power 
of  background  noise  is  about  30pWatts  with  a  maximum  sound  pressure  of 
3000pPascals  which  converted  into  decibel  units  gives  a  sound  pressure  level 
of  about  40dB.  However,  the  power  of  very  loud  music  may  be  30  Watts  with  a 
sound  pressure  of  3  Pascals  which  corresponds  to  a  sound  pressure  level  near 
lOOdB. 

The  transmission  loss  is  defined  as 


TL  =  -  10  log|  (170) 

where,  in  the  water  column,  I(r,z)  is  the  magnitude  of  the  acoustic  intensity 
in  Equation  (33)  and  is  the  reference  intensity  at  one  meter  from  the 
source  in  the  water  column.  Since  spherical  spreading  of  the  waves  occurs  at 
one  meter  from  the  source,  the  reference  intensity  is  equal  to  the  square  of 
the  time-averaged  rms  pressure  at  one  meter  from  the  source  divided  by  the 
acoustic  impedance  at  this  same  distance.  The  reference  pressure  is  given  by 
Equation  (28).  Substitution  of  the  scalar  potential  into  Equation  (28)  gives 

N  .. 

p(r,z)  =  -  w  p(z  )  p(z)  y  u  (z  )  u  (z)  H”’(k  r)  (171) 

4  0  LnOn  On 

n=  1 


56 


NSWC  TR  89-170 


and  this  expression  substituted  into  the  transmission  loss  expression  gives 


TL  (r,z)  =  20  log 

C 


(172) 


which  is  called  the  coherent  transmission  loss  because  the  phase  factor  of 
each  mode  has  been  taken  into  account  in  the  summation,  and  the  absolute  value 
is  taken  after  the  summation  is  completed.  The  layer  subscript  J  has  been 
erased  from  the  variables  to  minimize  the  complexity  of  the  equations.  Note 
that  the  transmission  loss  is  also  a  function  of  the  depth  of  the  source 
and  the  fact  that  the  acoustic  intensity  is  proportional  to  the  pressure 
squared  has  been  used.  This  coherent  transmission  loss  is  highly  variable  in 
space  due  to  phase-dependent  interference  effects  among  the  eigenfunctions  and 
a  smoother  function  is  more  appropriate  for  sonar  predictions. 

The  detailed  interference  effects  may  be  -veraged-out  to  yield  smooth 
transmission  loss  curves  by  summing  the  indiviaual  modal  energies.  The  result 
is  called  the  incoherent  transmission  loss.  The  resulting  incoherent 
transmission  loss  in  the  water  column  is 


TL  (r,z)  =  10  login^p^iz)  f  ju  (z  )  u  (z)  H''*(k  r)l^|  (173) 

1  I  '  n  0  n  0  n  1 

^  n=  1  ■' 

where  the  logarithm  is  of  base  10,  the  transmission  loss  is  always  purely 
real,  and  its  dimension  is  in  decibels. 

In  the  elastic  layers,  the  intensity  at  the  point  of  observation  is  now 
given  by  the  expression 


I(r,z)  =  I  P  .  ^  I  (174) 

where  P  is  the  stress  matrix  in  Equation  (111)  and  v  is  the  particle  velocity 
in  Equations  (55),  (60a),  and  (60b)  where  the  potentials  are  given  in 
Equations  (86)  and  (101).  With  these  equations,  the  vector  intensity  is 
written  as 
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t  =  P*v  =  I  r  +  I  z 


(175) 


where 


Hdw  /-V  5v  V  ,  ,  dw  9v  •.') 

.  X  (  ^  .  3^  )  ]  *  .  (  3^:  .  )j 


(176a) 


r  ,  V  5v  v  ,  f  dw  5v 

.=  |v.[(x*2m)  .  X  [  pi  »  3pi  ]  ]  .  ^  (  3j^  *  ]|  /i4,. 


(176b) 

Direct  substitution  of  the  scalar  potentials  in  Equations  (86)  and  (101)  into 
the  particle  velocity  components  in  Equations  (60a)  and  (60b)  gives 


V  (r,z)  =  y  fv  ) 

2  ^  '  2'  n 


(177) 


where  the  contribution  from  the  n  mode  is  given  by 


fv  1  =  (i/4)  p(z  )  u  (z  )  fu' (z)  +  V  (z)l  H*'’(k  r)  (178) 

'2^n  ono*-n  nn^o  n 


V  (r.z)  =  y  (v  ) 

r  x  r-'  n 


where 


(179) 


fv  1  =  (i/4)  p(z  )  u  (z  )  fu  (z)  +  v' (z)l  k  (k  r)  (180) 

'r'n  ono‘-n  n-'no  n 


and  the  prime  over  the  Hankel  function  stands  for  the  derivative  with  respect 
to  the  argument. 

By  the  same  token,  the  derivatives  of  the  components  of  the  particle 
velocity  are  given  by  the  summation  of  the  modal  contributions  given  by 
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(5v  /dr)  =  (t/4)  p(z  )  u  (z  )  fu' (z)  +  v  (z)]  k  H*^’'(k  r) 

z  n  ono^n  nn^no  n 


(181a) 


(5v  /Sr)  =  (t/4)  p(z  )  u  (z  )  [u  (z)  +  v' (z)]  k^  H*^^"(k  r) 
rn  ono*“n  n-*no  n 


(181b) 


(Sv  /Sz) 
z  n 


(i/4)  p(z  )  u  (z  ) 

o  no 


•(k^[u  (z)  +  v'  (z)]  -  k^(z)  u  (z)l  H 

^  n>-  n  n  J  n  J 


'"'(k  r) 

o  n 

(181c) 


(Sv  /Sz) 
r  n 


=  (i/4)  p(z  )  u  (z  )  fu' (z)  +  k^  V  (z)  -  K^{.z)  V  (z)''  k  H*'’'(k  r) 

Ono^*^  nn  nJnfi  n 


(181d) 


where  H*^*"(x)  =  -  H^^’(x)  -  H*^’'(x)  /  x.  After  calculating  the  particle 

O  0  0 

velocity  and  their  derivatives,  these  are  substituted  in  the  intensity  and 
this  one  in  the  transmission  loss  equation 


TL  (r.z)  =  5  log 

c  r  z 


(182) 


to  obtain  the  coherent  transmission  loss  in  the  elastic  bottom  layers  with 
unit  reference  intensity. 

For  the  incoherent  transmission  loss  in  the  solid  layers,  the  components 
of  the  intensity  are  calculated  for  each  mode,  i.e., 


(I  ) 

r  n 


(v  )  ,av 


[/ov  X  vv  ;  rUV  V 

^  [  6/  (^]j 


(183a) 
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V 


(v  ) 
r  n 


and  they  are  added  using 


/  iu 
(183b) 


and 


I 

r 


n=  1 


(184a) 


N 

I  =7(1) 

z  L,  z  n 
n  =  1 


(184b) 


to  substitute  into  Equations  (182),  where  the  subscript  in  the  transmission 
loss  is  replaced  by  the  incoherent  one. 

Use  caution  when  predicting  the  propagation  of  sound  using  any  incoherent 
transmission  loss  expression  since  this  is  just  an  approximation  to  obtain  a 
smooth  curve.  The  coherent  and  incoherent  transmission  loss  curve  should  be 
displayed  together  to  be  aware  of  the  variability  of  the  result. 

As  an  example,  the  transmission  loss  (see  Figure  10)  is  calculated  using 
the  seven  eigenvalues  in  Figure  8.  The  solid  curve  represents  the  coherent 
transmission  loss  and  the  dash  curve  is  the  incoherent  calculation.  The  depth 
of  the  source  and  the  receiver  is  112.0  meters.  The  high  oscillations  of  the 
coherent  curve  is  caused  by  the  constructive  and  destructive  interference  of 
the  trapped  modes.  The  radiating  modes  hardly  contribute  to  the  transmission 
loss  calculations  in  the  water  column  because  of  the  large  imaginary  part  of 
their  respective  eigenvalues. 

As  usual,  if  the  acoustic  properties  used  are  inaccurate,  then  the 
calculated  transmission  loss  will  be  erroneous. 

Hamilton  has  been  actively  involved  in  the  determination  of  the  geo¬ 
acoustic  properties  of  the  ocean  floor  by  the  use  of  more  easily  measurable 
quantities.  His  work  on  the  determination  of  the  compressional  sound  speed  in 
the  elastic  bottom  is  presented  in  References  66  and  67  while  some  results  of 
the  shear  wave  velocity  are  given  in  Reference  68.  The  shear  and  compressional 
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attenuation  coefficients  are  given  in  References  69  and  70  and  an  informative 
geo-acoustic  compendium  is  available  in  Reference  71.  Laboratory  measurements 

7  p  7  3 

of  the  shear  speed  and  the  shear  attenuation  coefficient  have  been  made  as 

functions  of  depth  and  frequency,  but  there  are  concerns  about  the 

effectiveness  of  a  laboratory  to  mimic  the  conditions  that  the  sediment 

encounters  under  the  unusual  pressure  and  temperature  of  an  ocean  column.  An 

accurate  non-destructive  method  must  be  created  for  the  proper  evaluation  of 

these  acoustic  properties.  Even  though  there  are  disagreements  about  most 

geo-acoustic  properties,  a  couple  of  inequalities  have  been  created  based  on 

2  8 

experimental  observations.  The  ratio  of  shear  to  compressional  attenuation 
coefficients  should  satisfy 


b  /c  s  V^TtS  (185a) 

J  J 

and  the  ratio  of  shear  to  compressional  attenuation  coefficients  should 
satisfy 

8^/a^  s  0.75  (185b) 

for  each  elastic . layer  from  j  =  J,  ...,  F-1. 
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CHAPTER  9 

COMPARISON  VIIR  SIMPLE  RANGE- INDEPENDENT  MODELS 


A  few  very  simple  range- independent  benchmark  ocean  models  will  be 
considered  to  compare  our  calculations. 

1.  The  first  model  is  a  one-layer  water  wave  guide  with  constant  acoustic 
properties,  a  pressure-release  surface,  and  a  rigid  bottom.  This  model 
resembles  the  infinite  well  in  quantum  mechanics  with  the  exception  that  the 
boundary  condition  at  the  bottom  is  u' (z=D)  =  0.  The  soft  surface  is  described 

n 

mathematically  by  u  (z=0)  =  0  and  the  eigenfunction  that  satisfies  the  wave 

n 

equation  and  both  boundary  conditions  is  given  by 


u  (z)  =  a  sin(7j  z),  n  =  1,  2,  3 . N  (186) 

n  n  n 


where 


7)  =  f  =  (n  -  l/2)7r/D  (187) 

which  provides  the  eigenvalues 

k^  =  <//c^  -  (n  -  1/2)V/D^  (188) 

n 

without  the  necessity  of  a  characteristic  equation  to  search  for  the  zeros. 
Note  that  the  fundamental  mode  is  given  by  n=l  and  as  n  increases,  the 
eigenvalue  decreases  towards  zero  until  a  limit  is  reached  when  k  converts 

n 

from  purely  real  (propagating  mode)  to  purely  imaginary.  The  maximum  number  of 
modes  N  is  obtained  by  setting  k  =0  giving 

N 
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N  =  —  +  i.  (189) 

nc  2 

The  normalization  equation  provides  the  amplitude  of  the  eigenfunctions.  The 
resulting  amplitude  is 


a^  =  ?./(pD)  (190) 

n 

where  p  is  the  density  of  the  water  and  the  subscript  n  is  not  necessary 
because  the  normalization  coefficient  is  constant  for  all  the  modes.  Finally, 
substitution  of  the  normalized  eigenfunctions  in  the  equation  for  the  coherent 
transmission  loss  gives 


TL(r.z)  =  20  log 


N 

-2ni/D  ^  sin[(n  -  1/2)7iz^/D]  sin[(n  -  1/2)tiz/D]  H^^^(k^r) 

n=  1 

(191) 


where  the  density  of  the  water  has  been  canceled  out  from  the  transmission 
loss  calculation  providing  no  contribution.  The  same  result  would  be  obtained 
if  the  orthonormalization  condition  has  no  weighting  function.  The  density 
would  contribute  only  if  a  finite  impedance  mismatch  exists  in  the  wave  guide. 
Note  also  that  this  model  allows  the  introduction  of  the  attenuation 
coefficient  as  the  imaginary  part  of  the  sound  speed.  This  provides  complex 
eigenvalues  but  the  eigenfunctions  are  still  purely  real. 

If  the  liquid  layer  is  200.0  meters  deep  and  has  a  sound  speed  of  1500.0 
m/s  then  the  25.0  hz  source  mentioned  in  the  previous  sections  would  excite 
seven  modes  with  the  eigenvalues  given  in  the  second  column  of  Table  2. 


2.  The  second  model  is  similar  to  the  first  one,  but  the  bottom  boundary  is 
a  pressure-release  interface.  In  this  case,  the  eigenfunction  is  still  given 
by  Equation  (186),  but  with 


=  (n  7i)  /  D  ,  n  =  1,  2,  3 . N  (192) 


where 
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N  =  (w  D)  /  (n  c)  (193) 

and  the  maximum  amplitude  of  the  eigenfunction  is  still  given  by  Equation 
(190).  The  eigenvalues  are  now  given  by 

1  2  2.2  2  2  ._2  f.nj,-, 

k  =  w  /c  -  n  n  /D  (194) 

n 

and  the  final  transmission  loss  equation  is 


TL(r.z)  =  20  log 


n 

-2nl/D  ^  sin[(nn2^)/D]  sin [(nrrz )/D]  H^^*(k  r) 


n  =  1 


.  (195) 


The  six  resulting  eigenvalues  for  the  case  of  a  200.0  meters  deep  water  column 
with  a  sound  speed  of  1500.0  m/s,  a  soft  bottom,  and  a  25.0  hertz  continuous 
wave  are  given  in  the  third  column  of  Table  2.  Note  that  the  soft-bottom 
eigenvalues  are  practically  located  half-way  between  the  locations  of  the 
rigid-bottom  eigenvalues  and  it  is  expected  that  the  true  eigenvalues  for  a 
wave  guide  with  penetrable  bottom  be  located  between  these  two  limiting  cases. 


3.  Another  way  to  compare  the  calculated  eigenvalues  is  to  use  the 
perturbation  method  for  small  attenuation  coefficients.  If  the  elastic  media 
has  a  negligible  shear  contribution  and  the  compressional  attenuation 
coefficient  in  each  layer  is  very  small,  then  both  methods  must  give  nearly 
the  same  answer. 

Consider  the  complex  eigenequation. 


—  u  (z)  +  [k^U)  -  u  (z)  =  0  (196) 

dz^  "  "  " 

where  we  redefine  the  wave  number  as  A(z)  =  k(z)  cea(z),  c  is  used  here  to 
keep  track  of  the  effects  of  every  term  in  the  resulting  approximate  complex 
eigenequation  and  it  will  be  set  to  unity  at  the  end  of  the  calculations,  a(z) 
is  the  attenuation  coefficient  in  nepers/meter,  and  k(z)  =  (j/c(z).  The 
complex  wave  number  in  Equation  (196)  makes  the  eigenvalues  and  eigenfunctions 
complex.  If  a(z)  <<  k(z),  then  we  can  use  the  perturbation  method  to  obtain  a 
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more  accurate  transmission  loss.  In  this  case  we  will  write 


,  ,  (0)  (1)  ^  2  (2) 
u  (z)  — »  u  +  c  u  +  e  u 

n  n  n  n 


(197a) 


.(0)  -d)  2.(2) 

X  +  e  X  ♦  E  X 


(197b) 


which  substituted  in  the  complex  eigenequation  gives, 


,2,  ,  2  2,  .  .(0)  -(1)  2.  (2)lf  (I 

-  +k  (2)+2u:k(2)a(z)-e  a  (z)-X  -eX  -e  X  u 

n  n  n  J  n 


(0)  (1)  2 

+  e  u  +  E 
n  n 


'S))  - 

u  =^0 

"  J 


(198) 


which  is  an  approximation  to  the  complex  eigenequation  Equation  (196)  due  to 
the  expansions  Equations  (197). 

Combining  the  c°  terms  of  this  equation  gives  the  order  solution  to 
the  problem,  or 


d  (0)  f,  2,  (0)  ,  (0)  „ 

-  u  Ik  (z)  -  X  ]  u  =  0 

,  2  n  n  n 

dz 


(199) 


which  is  the  unperturbed  eigenequation  that  has  been  solved  for  the  purely 

real  eigenvalues  x'°'  =  k^  and  eigenfunctions  u^°^  =  u  .  This  unperturbed 

n  n  n  n 

eigenequation  corresponds  to  Equation  (69). 

Combining  the  terms  with  c',  which  corresponds  to  the  first  order 
perturbation  terms,  gives 


d  (1)  r-,.,  ,  ^  r  \  ^  n  ^  <0)  ,  (1) 

-  u  +  l2tk(z)a(z)  -X  Ju  +[k(z)-X  ju  =0  (200) 

,  2  n  n  n  n  ii 

dz 


where  the  unperturbed  eigenfunctions  are  normalized  by 


r  ^  r  1  <o) ,  .  (0) ,  .  .  _ 

p(z)u  (z)u  (z)dz=6 
J  n  m  r 


(201) 


where  z^  is  the  depth  of  the  resilient  bottom  of  the  basement, 
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Multiplying  Equation  (200)  by  pu*°*  and 

n 


integrating  yields 


J 


E  ,2  (1) 

b  d  U 

(0)  n  ,  ^ 

pu  -  dz  + 

n  ,  2 

dz  0 


I 


pu  [2<.k(z)a(z)-X  ]u  dz  + 

n  n  n 


i 


(0)  ,,  2,  ,  V  (0)  1  (1)  ,  _ 

pu  Ik  (z)-A  ]u  dz=0 

n  n  n 


(202) 


where  using  the  orthonormality  condition  of  the  unperturbed  eigenfunctions  in 
the  second  term  of  this  equation,  integrating  by  parts  twice  the  first  term, 
and  using  the  boundary  conditions  at  every  interface  to  cancel  out  the  surface 
contributions  gives 


J 


t  .2  (0) 

b  d  U 

(1 )  n  .  T 

pu  -  dz  +  2t 

n  ,2 

dz  0 


pu  k(z)a(z)u  dz  + 

n  n 


i 


(0)  .,2,  ,  .  (0)  ,  (1)  ,  .  (1) 
pu  [k  (z)-X  ]u  dz  =  A 

n  n  n  n 


(203) 

and  with  the  help  of  Equation  (199)  the  first  and  second  integrals  cancel  out 
giving  us  the  expression 


A^^’=  2t  f  p  k(z)  a(z)  lu'°M^dz  (204) 

n  qJ  '  n  ' 

which  is  the  first  order  perturbation  term  for  the  eigenvalue  and  its  values 
are  purely  imaginary. 

Now  we  write  the  perturbed  part  of  the  eigenfunction  under  the  basis  of 
the  unperturbed  part  since  this  is  an  orthonormal  basis,  i.e.. 


(1)  r 


A  u 

run  m 


(0) 


(205) 


and  substitute  in  Equation  (200)  to  obtain 


,2  (0) 
d  u 


^  (0) 

(k  (z)-A 


[2ik(z)a(z)-A**^  ]u*°* 

n  n 


0 


(206) 


then  multiply  by  and  integrate  as  done  before.  Integration  by  parts 

twice  cancels  a  few  terms,  and  the  orthonormality  condition  yields 
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A  = 

nl 


2i 


,  (0) 

A  0 
1 


r 


,  ,  ^  ,  X  (0)  (0) 

p  k(z)a(z)  u  u  dz 

n  1 


(207) 


which  is  in  terms  of  the  unperturbed  eigenfvinctions  and  eigenvalues,  is 
directly  proportional  to  the  absorption  coefficient,  and  is  a  purely  imaginary 
term. 

In  the  cases  of  trapped  modes,  where  the  imaginary  part  of  the 

eigenvalues  is  extremely  small,  we  can  rely  on  the  rapid  convergence  of  the 

perturbation  method  and  forget  about  a  second  order  perturbation  term.  When 

radiating  modes  are  taken  into  account,  we  must  consider  calculating  the 

second  order  perturbation  term. 

2 

The  e  terms  of  Equation  (200)  into  a  second  order  equation  gives 


d  (2)^  .,2,  ,  (0)  1  (2)^  f'.-i  /■  ^  ^  ^  <o)  1  (1)  r  2,  ,  ^  (2)  ,  (0) 

-  u  +  [k  (z)-A  ]u  +  I2tk(z)a(z)-A  ]u  =  [a  (z)+A  ]u 

dz2  n  n  n  n  n  n  n 


which  multiplied  by  pu^°*  and  integrated  as  done  with  the  first  order 


(208) 


eigenvalue  leads  us  to  the  equation 

z  z  z 

-(2)  „.f  ,  (0)  (1).  ..(DT  2,  ,,  (0)|2, 

A  =  2t  pk(z)a(z)u  u  dz  -  A  pu  u  dz  -  pa  (z)  u  dz 

n  n  n  n  n  n  n 

(209) 

where  substituting  Equation  (205)  and  the  orthonormality  condition  of  the 
unperturbed  eigenfunctions  gives 


a'^'=  ZiY  A  rpk(z)a(z)u*°*u*°’dz  -  [pa 
n  Pi  nni  J  n  m  J 

Bt  0^ 


A  I 

(z)  u  dz 

'  n  ' 


(210) 


or  with  Equation  (207)  we  get  the  simpler  form 


y  (a'”’  a"”)  -  C 

n  P,  nm  n  m  J 

n'' 


2f  A!  (0)|2. 

(z )  u  dz 

'  n  ' 


(211) 


which  is  purely  real  and  a  much  smaller  term  since  it  is  proportional  to  a  . 
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If  we  write 

(2)  ^ 

u  =  )  B  u 

n  nm  m 

m 

then  Equation  (208)  becomes 


(212) 


d  (0)X  _  f.2  .(0),  (0)  r  ,  .  (1),  (0)  r  2  ^(2),  (0) 

B  -  u  +)  B  [k  -A  Ju  +)  A  [2tka-A  Ju  =  [a  -A  Ju 

nis,2mLinffl  nm^nm  nm  nn 

m  QZ  m  ID 

(213) 

which  multiplied  by  and  integrated  using  integration  by  parts  and  the 

orthonormality  condition  reduces  the  equation  to 


(a‘°’-  a‘°’)  B  =  a'^’a  -  y  a  a  .  l*n  (214) 

1  n  nl  n  nl  ^  nxn  Im 

m 


which  makes  B  purely  real  and  directly  proportional  to  a^. 

nl 

We  have  already  assumed  layers  of  constant  density  in  order  to  simplify 
the  elastic  wave  equation.  Therefore,  we  may  define  an  element  of  a  G-matrix 
as 

J  + 1  Z 

I  '  '  -  j  1 

G  =  2i  )  p  I  k  (z)  a  (z)  u'°’(z)  u'°’(z)  dz  =  G  (215a) 

nm  / _ ^  J  J  J  J  nj  mj  mn 

J  =  1  J 

and  that  of  an  H-vector  as 


1=1 


(215b) 


Note  that  all  elements  of  the  G-matrix  are  purely  imaginary  and  symmetric, 
while  those  of  the  H-vector  are  r..rely  real.  These  integrals  must  be 
evaluated  in  order  to  calculate  the  perturbed  parts  of  the  eigenvalues  and 
eigenf  unctions. 

Now  the  first  order  perturbation  term  of  the  eigenvalue,  Equation  (204), 
becomes 

a‘^’=  G  (216) 

n  nn 
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which  tells  us  that  the  diagonal  components  of  the  G-matrix  are  the 
first-order  perturbation  term  of  the  eigenvalues.  The  second  order  term, 
Equation  (209),  simplifies  to 


X 


(2) 


n 


r 

l*n 


^(0)_  ^(0, 

1 


n 


(217) 


which  substituted  into 


c.  x“’+  x'^’  (218) 

n  n  n  n 

gives  the  perturbed  eigenvalues  of  the  problem.  Since  is  the  only 

n 

contributor  to  the  imaginary  part  of  the  eigenvalue,  we  may  define 


3?^  =  X‘°^  X^^’ 

n  n  n 


(219) 


as  the  real  part  of  the  eigenvalue.  Then  to  obtain  k  from  Equation  (218)  we 

n 

expand  its  square  root  as  follows: 


k  ^ 

n 


( :  ) 


2 


( 1  ) 


(220) 


2  (1 ) 

where  we  have  assumed  that  3?  >>  X  .  Now  the  imaginary  part  of  the 

n  n 

eigenvalue  will  be  defined  as 


3  ^ 

n 


n 


2  3? 


n 


(221) 


which  is  the  same  expression  in  Equation  (474)  of  Reference  53  where  only 
first  order  perturbation  has  been  used.  By  the  same  token,  the  real  part  is 
given  from  Equation  (219) 
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A' 

3?  c.  -  (222) 

n  n  2v/x'°’ 
n 

where  it  is  assumed  that  A^^*  <<  A*°'  and  the  same  power  expansion  has  been 

n  n 

used.  Equation  (221)  is  a  crude  approximation  made  by  many  underwater 
acousticians  and  it  can  be  avoided  by  taking  the  complex  square  root  of 
Equation  (218). 

As  the  first  order  correction  of  the  eigenfunction,  Equation  (207)  simply 
becomes 


A  = 


n  1 


nl  .(0)  .(0) 

A  A 
n  1 


,  n^l 


(223) 


and  for  the  second  order  correction,  Equation  (214),  we  get 


B  = 

nl 


-  G  G 

nn  nl 

(a‘°>-  a‘°>)^ 

n  1 


G  G 

njn  tn  1 


a‘'^>)(a'°’-  a‘°’) 

m*n  n  m  n  1 


n»l  (224) 


which  substituted  into  the  equation 


,  ,  (0)  \  ,  (0) 
u  (z)  =  u  +  )  [A+B  ]u 

n  n  /  nl  n  1  1 


(225) 


1  *n 


gives  a  better  estimate  of  the  eigenfunction. 

It  is  left  to  properly  evaluate  the  G-matrix  and  the  H-vector  in 
Equations  (215).  They  can  be  obtained  by  numerical  integration  or  by  the 
approximate  method  developed  in  Reference  53.  Since  we  are  interested  here  in 
the  simple  case  of  a  semi-infinite  fluid-type  bottom  Equations  (215)  become 
trivial  integrations  of  exponential  functions.  The  resulting  trapped 
eigenvalues  are  given  in  the  fourth  column  of  Table  2.  Note  that  the  second 
order  perturbation  term  was  not  enough  to  make  the  real  part  of  the  eigenvalue 

closer  to  the  exact  one.  Therefore,  the  attenuation  coefficient  chosen  by 

4  2  4  3 

Miller  ’  is  too  high  for  the  perturbation  method  to  produce  accurate 
results  for  the  trapped  eigenvalues,  and  it  is  presumably  worst  for  the 
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radiating  ones. 

4.  The  next  model  is  a  layer  of  water  over  a  semi-infinite  elastic  layer 
which  was  used  by  Ellis^^  to  model  the  propagation  of  underwater  explosives 
over  an  ocean  floor  made  of  chalk.  Both  layers  have  constant  acoustic 
properties.  Since  the  surface  of  the  water  is  pressure-release,  the 
eigenfunction  is  given  by 

u  (z)  =  a  sin(T)  z)  (226) 

nl  1  nl 

2  2  2 

where  t)  =  k  -  k  and  where  z  =  0  i  z  <  z  .  The  compress ional  and  shear 
eigenfunctions  in  the  semi-infinite  layer  radiated  without  reflecting  back, 
therefore  the  compressional  mode  is  given  by 


u  (z)  =  a^  expilri  z)  (227) 

n2  2  n2 

2  2  2 

where  t)  =  k  -  k  and  the  shear  mode  is 
n2  2  n 


V  (z)  =  b  expCiy  z)  (228) 

n  n 

2  2  2 

where  y  =  k  -  k  .  The  three  unknown  constants,  a  ,  a  ,  and  b,  are  determined 
n  n  12 

by  the  three  liquid-solid  boundary  conditions  in  Equations  (113),  (115),  and 
(116).  Direct  substitution  of  the  given  eigenfunctions  in  the  three  boundary 
conditions  provides  two  equations  with  the  three  unknowns,  i.e., 


2 

a  T?  cos(7}  z  )  =  a  IT)  expdr?  z  )  +  k  b  expdy  z)  (229a) 

1  nl  nl  2  2  n2  n2  2  n  ^  n2 

and 

2  2  2  2 
p  K  a  sin(T)  z  )  =  p  (k  -  2k  )a  expdij  z  )  -  2p  k  h  iy  expiiy  z  ) 

1  1  nl  2  2  n  2  n2  2  2  n  n2  n2  2 

(229b) 

and  a  third  equation  with  two  of  the  unknowns,  i.e.. 


2  a  IT)  exp(<T/  z  )  +  (2k^  -  k^)  b  exp(<>  z  )  =  0. 
2  n2  n2  2  n  ^  n2  2 
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If  this  third  equation  is  solved  for  one  of  the  unknowns  to  be  substituted 
into  the  other  two  equations,  we  obtain  two  equations  with  two  unknowns.  The 
two  equations  are  divided  to  eliminate  the  remaining  unknowns,  and  to  form  the 
characteristic  equation 

W(k  )  =  7)  tan(T}  z  )  +  i  p  /p  f  (2k^  -  +  4y  7)  k^']=0 

n  n2  nl  2  ^^2  1  >-  n  n2  n2  n  -I 

(230) 

which  is  a  complex  equation  even  if  the  attenuation  coefficients  are  not 
included  in  the  model.  The  transmission  loss  is  calculated  using  the  same 
equations  derived  in  the  previous  section. 

Substitution  of  the  parameters  in  the  "Fluid  bottom"  case,  shown  in  Table 
1,  into  this  simple  model  gives  the  three  eigenvalues  displayed  in  the  last 
column  of  Table  2.  These  eigenvalues  agree  with  double  precision  accuracy  with 
the  ones  obtained  using  our  multilayered  model,  hence  the  figures  in  the  last 
column  of  Table  2  represent  the  solutions  from  both  methods.  Double  precision 
accuracy  was  also  found  between  the  resulting  eigenvalues  from  this  simple 
model  and  our  multilayer  model  when  the  bottom  is  a  semi-infinite  layer  of 
Clay-silt,  Sand,  Basalt,  or  Chalk  (see  Table  1). 

Note  that  these  eigenvalues  are  located  between  the  rigid-bottom  and 
soft-bottom  eigenvalues.  Hence,  these  two  simple  cases  may  be  used  to 
establish  limits  to  the  eigenvalues  searched. 

With  the  use  of  these  simple  models,  the  range-independent  multilayer 
model  has  been  compared  in  the  limit  when  the  number  of  layers  is  a  minimum 
(with  or  without  shear  waves).  Therefore,  some  transmission  loss  results  will 
be  made  for  the  various  bottom  types  given  in  Table  1. 

The  first  case  to  consider  is  a  water  column  2C0.0  meters  deep  over  the 
semi-infinite  sand  ocean  floor  with  the  acoustic  properties  in  Table  1.  A 
total  of  seven  eigenvalues  were  found  (see  Figure  11)  for  a  25.0  Hz  source. 

The  real  part  of  the  normalized  fundamental  mode  is  shown  in  Figure  12  where 
the  solid  curve  represents  the  compressional  eigenfunction  and  the  dotted 
curve  is  the  shear  contribution.  The  imaginary  part  of  this  fundamental  mode 
is  displayed  in  Figure  13  where  the  discontinuity  at  the  liquid-solid 
interface  is  caused  by  the  impedance  mismatch  and  the  presence  of  the  shear 
contribution.  The  first  three  eigenfunctions  are  T-R  modes  because  the 
compressional  sound  speed  in  the  basement  layer  is  greater  than  the  sound 
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speed  in  the  water  column  (1500.0  m/s)  and  the  shear  sound  speed  is  smaller 
than  the  water  sound  speed.  The  fourth  and  higher  depth  functions  are  R-R 
modes.  The  real  part  of  the  fourth  mode  is  given  in  Figure  14  where  the 
compressional  eigenfunction  starts  to  oscillate  into  the  bottom.  The  imaginary 
part  of  this  mode  is  plotted  in  Figure  15.  Note  that  the  first  and  fourth  mode 
are  damped  out  when  they  reach  the  depth  of  600.0  meters.  At  this  depth  the 
higher  order  radiating  modes  take  over  in  the  transmission  loss  calculation. 

For  best  visualization  of  the  propagation  of  sound  in  range  and  depth,  the 
three-dimensional  plots  and  the  contour  plots  have  been  proven  to  be  useful 
tools.  The  advantage  of  the  three  dimensional  plots  is  that  every  calculated 
point  is  plotted.  However,  the  disadvantage  is  that  it  is  more  difficult  to 
visually  extrapolate  the  numerical  values  of  any  point.  The  advantage  of  the 
contour  plot  is  that  it  is  a  two-dimensional  plot  and  numerical  values  can  be 
roughly  extrapolated  visually,  but  not  all  the  calculated  values  are  used  to 
obtain  the  contours  and  there  is  less  information  in  this  type  of  graphical 
display.  For  the  various  tastes  of  the  readers,  both  plots  will  be  displayed. 

The  three-dimensional  transmission  loss  plot  versus  range  (in  kilometers) 
and  depth  (in  meters)  is  given  in  Figure  16  for  the  same  case  of  the 
semi-infinite  sand  basement  and  a  25.0  hertz  source  located  100.0  meters  deep. 
The  transmission  loss  at  the  liquid-solid  interface  has  been  intentionally 
omitted  from  the  plot  to  mark  the  location  of  the  water  depth.  Note  the 
oscillatory  behavior  of  the  transmission  loss  surface  at  short  ranges  which  is 
caused  by  the  interference  of  the  seven  normal  modes  of  this  wave  guide. 
However,  the  high  imaginary  component  of  the  excited  eigenvalues  (see  Figure 
11)  causes  all  the  modes,  except  the  fundamental  mode,  to  completely  damp  out 
at  ranges  beyond  10.0  kilometers.  The  contour  plot  for  the  results  in  Figure 
16  is  displayed  in  Figure  17.  The  legend  at  the  bottom  of  the  contour  plot 
relates  the  curve  type  to  the  contour  transmission  loss  value  in  decibels. 

The  next  case  to  consider  is  the  200.0  meters  deep  water  column  over  a 
semi-infinite  clay-silt  ocean  floor  with  the  acoustic  properties  in  Table  1.  A 
total  of  seven  eigenvalues  were  found  (see  Figure  18)  for  the  25.0  Hz  source. 
The  real  (imaginary)  part  of  the  calculated  fundamental  mode  is  plotted  in 
Figure  19  (see  Figure  20  also).  A  very  interesting  observation  made  is  that 
all  the  modes  are  R-R  modes  which  radiate  into  the  bottom.  This  is  because  the 
compressional  speed  in  the  solid  basement  is  almost  equal  to  the  water  sound 
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speed.  Also,  there  is  hardly  any  contribution  from  the  shear  waves  because  of 
the  relatively  low  shear  speed  of  the  clay-silt  ocean  floor.  Since  the  excited 
modes  display  the  same  oscillatory  behavior,  they  are  not  illustrated.  The 
three-dimensional  transmission  loss  plot  in  Figure  21  exhibit  a  faster  decay 
of  the  contribution  from  the  excited  modes,  compared  to  the  sand  bottom  case 
in  Figure  16,  despite  the  lower  attenuation  coefficients.  At  ranges  greater 
than  5.0  kilometers  the  fundamental  mode  becomes  the  only  contributor  to  the 
transmission  loss.  This  rapid  decay  of  the  higher  modes  is  caused  by  the  very 
low  compressional  sound  speed  of  the  clay-silt  basement.  Only  R-R  modes  are 
excited  by  the  source  and  these  radiating  modes  are  strongly  affected  by  the 
acoustic  absorption  of  the  bottom.  The  trapped  modes  in  the  sand  bottom  case 
are  evanescent  in  the  bottom  and  are  much  less  affected  by  the  high 
attenuation  of  sand.  Comparison  of  the  contour  plot  for  clay-silt  bottom  (see 
Figure  22)  with  the  one  for  sand  bottom  (see  Figure  17)  shows  that  the 
fundamental  mode  propagates  further  in  the  water  column  if  the  bottom  is  made 
of  sand.  The  same  conclusion  can  be  obtained  by  comparison  of  Figures  16  and 
21,  but  it  is  easier  to  visualize  with  the  contour  plot. 

It  is  almost  impossible  to  find  a  water  column  directly  over  basalt. 
However,  the  case  will  be  considered  only  for  its  interesting  acoustic 
features.  Note,  from  Table  1,  that  basalt  has  lower  attenuation  coefficients 
relative  to  those  of  sand  or  clay-silt  and  that  both  sound  speeds  are  much 
higher  than  the  water  sound  speed.  The  eight  complex  eigenvalues  found  are 
plotted  in  Figure  23.  The  first  five  modes  are  T-T  modes  because  of  the  high 
sound  speeds  of  basalt.  Figure  24  displays  the  real  part  of  the  normalized 
fundamental  depth  function  and  Figure  25  is  the  fifth  normal  mode?  As  the  mode 
number  increases,  the  effects  of  bottom  absorption  increases  and  the  imaginary 
part  of  the  eigenvalue  in  Figure  23  increases.  The  sixth  mode  (see  Figure  26) 
is  the  first  R-T  mode  and  absorptive  effect  to  the  oscillatory  compressional 
mode  is  different  from  the  effect  to  the  exponential  one,  hence  the  irregular 
pattern  displayed  in  Figure  23.  The  seventh  eigenvalue  corresponds  to  the 
second  R-T  mode  (see  Figure  27),  and  the  last  eigenvalue  represents  the  only 


■•^To  minimize  the  number  of  figures  displayed  for  this  case,  the  imaginary  part 
of  the  eigenfunctions  are  omitted. 
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R-R  mode  (see  Figure  28)  of  the  set. 

The  three-dimensional  plot  of  the  coherent  transmission  loss  for  this 
case  is  given  in  Figure  29.  All  the  modes  contribute  to  the  transmission  loss 
for  ranges  smaller  than  5.0  kilometers.  In  the  case  of  larger  ranges,  only  the 
first  four  normal  modes  are  needed.  The  second,  third,  and  fourth  modes  are 
not  damped  out  at  10.0  kilometers  because  they  are  T-T  modes  that  propagate 
mostly  in  the  water  column  and  experience  a  negligible  effect  from  the 
absorption  of  the  bottom.  The  fifth  mode  is  also  a  T-T  mode,  but  the  imaginary 
component  of  its  eigenvalue  (see  Figure  23)  is  much  higher,  causing  the  mode 
to  dampen  out  rapidly. 

The  contour  plot  for  this  same  data  is  provided  in  Figure  30  where  the 
discontinuity  in  the  transmission  loss  at  the  liquid-solid  interface  is  a 
result  of  the  impedance  mismatch  of  the  boundary.  The  highly  oscillatory 
behavior  of  the  transmission  loss  makes  the  contour  plot  somewhat  complicated 
to  interpret.  In  this  case,  the  three-dimensional  plot  may  be  more  useful. 
However,  the  contour  plot  does  show  the  bundle  of  acoustic  energy  that 
scatters  the  bottom  several  times  at  the  critical  angle  of  incidence.  There 
are  17  surface  and  bottom  bounces  for  ranges  between  3.83  kilometers  and  14.85 
kilometers.  This  corresponds  to  a  critical  angle  of  incidence  of  about  72.8 
degrees  relative  to  the  vertical  axis.  This  plot  provides  a  relationship 
between  the  normal  mode  theory  and  the  ray  theory. 
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CHAPTER  10 

ADIABATIC  NORMAL-MODE  THEORY  WITH  SHEAR  HAVES 


After  verifying  that  the  range- independent  normal  mode  calculations  agree 
with  the  calculation  obtained  with  simple  benchmark  models,  the  next  step  is 
to  include  range  dependence  with  the  shear  effects  of  the  elastic  bottom 
sediments. 

The  range-dependent  Helmholtz  equation  in  the  water  layers  of  the  ocean 
wave  guide  is  given  as, 

+  k^(r,2)  v(r,z)  =  6(r)  6(z-z^)  (231) 

2 

where  now  k  displays  the  range-dependence  of  the  acoustic  properties  of  the 
ocean.  The  range-dependence  of  the  boundaries  are  displayed  in  the  boundary 
conditions  themselves. 

The  range-independent  solution  was  found  to  be  given  by, 

N 

<p{r,z)  =  -  p{z  )  y  u  (z  )  u  (z)  H'^’(k  r).  (232) 

4  oLnOn  0  n 

n=  1 

However,  in  the  range-dependent  case,  the  eigenfunctions  and  eigenvalues  vary 
with  range,  therefore  the  solution  may  be  written  using  the  quasi-  separation 
of  variables  as. 


N 

^(r,z)Hy  f(r)u(r,z)  (233) 

L  r\  n 

n  =  1 

where  u  (r.z)  are  taken  as  the  basis  depth  functions  that  satisfy  the  equation 
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- —  u  (r,z)  +  fk^Cr.z)  -  k^{r)]  u  (r,z)  =  0  (234) 

n  L  n  J  n 

and  the  orthonormality  condition 


z)  u  (r,z)  u  (r,z)  dz  = 

n  D 


5 

nin 


Direct  substitution  of  Equation  (232)  into  Equation  (231)  gives, 


(235) 


■|7^[f^(r)  u^(r,z)]  +  f^(r)[aVaz^  +  k^(r,  z)]u^(r,  z)j-  =  5(r)  5(z- 


^0^ 

(236) 


where 


=  i  ^  fr 

r  r  ar  ar 

V  J 


(237) 


and  we  may  substitute 


V^[f  (r)  u  (r,z)] 
r  ^  n  n 


u  V^f  +  2^f-^u  +  f  V^u 

nrn  rn  rn  nrn 


(238) 


to  obtain 


].{f 


(V  +  k  )f  lu 

r  n  n*^  r 


2  ^  f  •  ^  u 


+  f  V^u 

r  n  nrn 


-1 

2Trr 


5(r)  6(z-z 


(239) 


which  multiplied  both  sides  by  p(z)  u  (O.z)  and  integrated  throughout  depth 

n 

gives  the  inhomogeneous  range  equation 


where 


+  k^(r)lf 

n  n 


(r)  =  = —  6(r)  p(z  )  u  (O.z  ] 
2nr  On  0 


n  n 

2  y  f' (r)U  (r )  -  y  f  (r)W  (r) 

Lm  m  nm  Li  m  nrn 


fn=  I 


m  =  1 


(240) 


77 


NSWC  TR  89-170 


1)  (r)  £  r  ^p(z)  u  (r.z)  ^  u  (r.z)  dz  (241a) 

nm  J  m  OT  n 

and 

W  (r)  s  r  ^p(z)  u  (r,2)  -  I-  [r  u  (r,z)  dz  (241b) 

nm  m  r  or  Or  J  n 

are  the  elements  of  the  coupling  matrices  that  will  take  care  of  the  exchange 
of  energy  of  the  normal  modes  in  the  range-dependent  environment.  The  prime 
stands  for  the  derivative  with  respect  to  the  argument. 

In  the  case  where  the  acoustic  properties  and  the  boundaries  of  the  ocean 
wave  guide  slowly  vary  with  range,  the  coupling  integrals  are  negligible  and 
the  adiabatic  approximation  is  feasible.  The  adiabatic  range  equation  is 

r  +  k^(r)]f  (r)  =  6(r)  p(z  )  u  (O.z  )  (242) 

where  the  range-dependent  waveguide  will  be  divided  into  M  number  of  range- 
independent  segments.  The  procedure  is  to  calculate  a  fixed  number  of  trapped 
and  radiating  modes  for  each  range-independent  segment.  The  resulting  set  of 
eigenfunctions  provides  the  function  u  (r.z).  The  unknown  function,  f  (r),  is 

n  n 

obtained  from  the  range  equation  and  the  range  boundary  conditions. 

Range  segment  #1  is  defined  as  the  one  where  the  source  is  located.  The 
homogeneous  solution  for  the  first  range  segment  is 

V  (r)  £  H"’(^k  r)  +  *8  H‘^’(^k  r)  (243) 

n  n  0  n  n  0  n 

where  the  left-side  superscripts  is  the  range  segment  number,  and  the  unknown 
constants  are  to  be  determined  with  the  range  boundary  conditions. 

In  the  limit  as  ^k  r  — >  0  the  asymptotic  forms  of  the  Hankel  functions 

n 

H^'’(^k  r)  =  -  H'^’('k  r)  — »  Zi/n  log  ('k  r)  (244) 

On  On  e  n 

gives 

^f  (r)  Zi/n  Oot  -  V  )  log  {'k  r).  (245) 

n  n  n  e  n 

The  particular  solution  is  obtained  by  integrating  the  inhomogeneous  range 
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equation  over  a  small  cylinder  of  radius  a  containing  the  source,  i.e., 


.a  .a  .  .a 

I  ^f"(r)  dr  +  r  -  V' (r)  dr  +  T  ^k^(r)  (r)  dr  = 

qJ  "  qJ  r  n  O'*  " 

-p(z  )  u  (0,z  )/(2n)  r  dr. 

0  n  0  Jr 

0 

Integrating  by  parts  gives, 


(246) 


^f^(r)|°  +  ^f^(r)/r|^  ^  I  [  ~ 

-p(z  )  u  (0,z  )/(2n)  r  dr 

^  0  n  0  J  r 

r\*^ 


(247) 


which  in  the  limit  as  a  — >  0,  only  the  slope  at  r  =  0  and  the  integral  over 
the  delta  function  remains,  i.e., 


d'f  (r)/dr  —>  -  p(z  )  u  (0,z  )  /  (27ir) 

n  0  n  0 


which  yields 


V  (r)  — )  ^  p(z  )  u  (0,z  )  log  (^)c  r), 

n  271  0  n  0  e  n 


Equating  the  two  solutions  provides  the  relationship 


-  ^8  =  j  p(z  )  u  (0,z  ) 

n  n  4  0  n  0 


(248) 


(249) 


(250) 


where  the  right-hand-side  term  is  the  constant  in  the  range-independent 
solution  and  this  equation  will  be  used  as  the  relationship  between  both 
unknowns  in  the  first  range  segment. 

The  range  segments  labeled  2  to  M-1  are  characterized  by  the  homogeneous 
range  equation,  therefore  the  solutions  are. 


n  M  /  \  in  f « ( 1 X  D)  ^  f  •  ( 2 )  /  ^  I  \ 

f  (r)  =  a  H  (  k  r)  +  8  H  (  k  r) 

n  n  0  n  n  0  n 


(251) 
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where  2  s  m  <  M. 

The  last  range  segment  is  assumed  semi-infinite  and  with  no  source 
present,  therefore  only  the  divergent  solution  of  the  homogeneous  range 
equation  satisfies  causality.  The  solution  for  the  range- independent 
segment  is 


f(r)=  a  H  (kr). 

n  n  0  n 

The  unknowns  are  determined  by  the  radial  boundary  conditions. 


(252) 

These  are: 


1.  Continuity  of  the  normal  particle  velocity 

N  ^  N 


n=  1 


^  y  “f  (r)  u  (r,z)|  =  f;:  y  u  (r,2)| 

or  L  n  n  •  r=r  OT  L  n  n  '  r=r 

m  n=  1 


(253a) 


2.  Continuity  of  the  pressure 

N 


py'"f(r)u(r,z)  =py  ^  (r,  z)  I 

U  r\  n  r=r  C  n  n  '  r=r 


(253b) 


n  =  l  m  n=  1  m 

Since  the  eigenfunctions  already  satisfy  the  boundary  conditions  at  every 
range  and  depth,  in  a  slowly  varying  range-dependent  environment,  the 
conditions  to  satisfy  for  each  mode  are 


""f  (r)  =  continuous  (254a) 

n 

and 

'"f'(r)  =  continuous  (254b) 

n 

for  1  s  m  <  M. 

Application  of  these  radial  boundary  conditions  to  the  M-1  interface 

gives 


a  H  (  k  r  ) 

n  0  n  M-l 


.  H‘^>("-^k  r  )  = 

n  0  n  M-l 


a  H  (  k  r  ) 

n  0  n  M-l 


(255a) 


and 
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H- 


k  aH  (  kr  )+  PH  (  kr  )  =akH  (kr  ) 

n  n  1  n  M-1  n  1  n  M-1  J  n  n  1  n  M-1 


(255b) 


which  is  rewritten  in  matrix  form  as 


where 


and 


M-1 


W 


M-1 , 


n,M-l 


=  *'lH 


n,M 


M, 

A 

•1  n 


(256) 


'W  £ 

n,  1 


H‘‘’(^k  r  ) 
0  n  1 


H‘^’(-’k  r  ) 

0  n  i 


Jk  H”>('k  r  )  Jk  H'^^Jk  r  ) 

n  1  n  1  n  1  n  i 


n 


(257) 


(258) 


for  i,  1  =  2,  3,  4,  .  .  .  ,  M  and  where  =  0.  To  obtain  ^A  in  terms  of 

n  n 

A  ,  we  write 

n 


n 


,M-1 


(H 

n,M-l 


M 


n,M-l 


(259) 


The  M-2  boundary  has  the  relationship 


n,M-2  n  n,M-2  n 


(260) 


and  substituting  the  previous  relationship  for 


gives 


*'■2^  =  ('’■‘iH  )■'  "iH  "a  .  (261) 

n  n,M-2  n,M-2  n,M-l  n,M-l  n 

By  the  same  token,  we  can  propagate  the  solution  to  the  first  segment  with 


’iH  ‘a  =  ^IH  ^A  (262) 

n,  1  n  n,  1  n 


which  gives 
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V  (“w  )-' 

"  n,m  n,B 

ID  =  1 


(263) 


where  the  term  in  the  parenthesis  is  now  defined  as  the  2x2  matrix  X  and  the 
coefficients  in  the  first  range  segment  are  given  by 


and 


V  =  X  **0  (264) 

n  11  n 

=  X  "a  .  (265) 

n  21  n 


Substitution  of  these  relationships  into  Equation  (250)  gives 


u  (0, z  )  =  (X 

n  0  11 


21  ' 


H 

a 

n 


(266) 


which  is  solved  for  the  unknown  at  the  semi-infinite  range  segment  and  this 
solution  can  be  propagated  to  obtain  the  other  unknowns. 

The  potentials  for  the  range-independent  solid  layers  are  given  bv. 


and 


^(r.z)  =  [  [  ^  P(z^)  u^(z^) 

n=  1 

0(r.z)  =l[{  P(z^)  u^(z^) 

n=  1 


h”*  (k  r  )1  u  (z) 

0  n  n 

h''*  (k  r  ll  V  (z) 
0  n  J  n 


(267a) 


(267b) 


where  the  term  in  the  brackets  is  common  to  both  solutions.  With  this 
observation  in  mind,  the  range-dependent  solutions  in  the  solid  layers  will  be 
written  as 


and 


<pir,z) 


u  (r , z) 


(268a) 


i/»(r,z) 


N 


(r) 


V  (r,z) 

n 


(268b) 
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which  must  satisfy  the  four  boundary  conditions  throughout  range  and  depth. 

The  range  function  for  these  potentials  is  the  same  as  the  one  given  for 
each  range  segment  in  the  liquid  layers.  The  unknown  constants  are  evaluated 
by  the  four  radial  boundary  conditions.  One  of  them  is  that  the  tangential 
component  of  the  particle  velocity  must  be  continuous,  i.e., 

v(r,2)  =  m  =  continuous  (269) 

z  m  Id*:  n  ir 

^  '  m 

where  substitution  of  <f  and  0  gives 

f  (r  )  fu' (r  ,z}  +  V  (r  ,z)]  =  continuous.  (270) 

nm'-nm  nnm 

However,  since  the  eigenfunctions  already  satisfy  the  boundary  conditions  in 
Equation  (119),  then  all  we  have  left  to  satisfy  is  the  continuity  of  f  (r  ). 

n  ifi 

The  normal  component  of  the  particle  velocity  is  another  boundary 
condition  to  be  satisfied.  This  is  given  by 

V  (r  ,z)  =  [^(r  ,z)  +  ^  0(r  ,z)l  =  continuous  (271) 

r  m  dr  I  m  dZ  m  J 

and  substitution  of  ^  and  t}j  gives 

f'  (r  )  fu  +  v'l  +  f  (r  )  ^  fu  +  v'1  =  continuous  (272) 

nm^n  nin  OV  ^  n  T 

m  m 

where,  in  a  slowly  varying  environment,  the  change  in  (u  +v'  )  with  respect  to 

n  n 

range  is  much  smaller  than  the  change  of  the  Hankel  functions  in  f  (r)  with 

n 

range.  Therefore,  the  only  functions  to  make  continuous  are  f  (r)  and  f'(r). 

n  n 

The  same  conditions  are  found  from  the  continuity  of  P  and  P  .  With  this 

rz  zz 

adiabatic  approximation,  the  need  to  match  four  boundary  conditions  explicitly 
has  been  avoided  and  only  two  equations  must  be  satisfied.  The  equations  to 
match  turn  out  to  be  the  same  as  the  ones  in  the  liquid  layers,  therefore  the 
same  function  f  (r)  can  be  used  for  both  states  of  the  matter.  This  property 

n 

may  decrease  the  computation  time  by  orders  of  magnitude. 

After  all  the  coefficients  of  the  range-dependent  waveguide  are 
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determined,  the  coherent  and  incoherent  transmission  loss  in  the  solid  layers 
are  obtained  by  Equations  (176)  through  (184)  where  the  components  of  the 
particle  velocity  are 


and 


V  =  y  f  (r)  ru'(r,z)  +  V  (r,z)l  (273a) 

zt-n  ‘-n  nn  J 

n  =  1 

N 

V  y  f' (r)  fu  (r,z)  +  v'(r,z)'|  (273b) 

r  Li  n  n  n 

n  =  1 


where  it  has  been  assumed  that  the  change  of  the  eigenfunctions  with  range  is 
negligible  compared  to  the  change  in  the  Hankel  functions.  The  derivatives 
with  respect  to  range  are  given  by 


V  =  y  f'(r)  ru'(r,z)  +  V  (r,z)]  (274a) 

or  zLn  '-n  nn 

n  =  1 

and 

P  N 

—  V  y  f"(r)  fu  (r,z)  +  v'(r,z)l  (274b) 

dr  r  L  n  •-n  n  J 

n=l 

where  the  homogeneous  range  equation  gives 


f"(r)  =  -  f'(r)/r  -  k^(r)  f  (r).  (275) 

n  n  n  n 

The  derivatives  with  respect  to  depth  are  given  by 


and 


N 


n=  1 


(u  +v'  ) 
n  n 


(276a) 


“I 


n=  1 


f'  (r) 

n 


fu'  +  V  1 

^  n  n  n-' 


(276b) 


The  coherent  transmission  loss  in  the  range-dependent  water  column  is 
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TL  (r,z)  =  - 

C 


=  -20  log[  4,r  ] 


(277) 


and  the  incoherent  transmission  loss  is 


2  N 

TL,(r.z)  .  -10  log[(4.t  e{|^]  ] 


(278) 


where  the  range  function  f  (r)  and  the  eigenfunctions  are  complex  and  the 

n 

transmission  loss  is  real. 
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CHAPTET!  11 

COMPARISON  WITH  EXPERIMENTAL  MEASUREMENTS 


The  solutions  from  this  range-dependent  model  can  be  compared  with  other 
range -dependent  models  with  the  purpose  of  validating  its  results.  However, 
there  is  no  other  range-dependent  model  that  can  include  the  effects  of  shear 
waves  from  the  ocean  bottom.  Anyway,  if  there  were  an  opportunity  for  the 
inter-model  comparison,  their  agreement  Coes  not  rule  out  the  possibility  tnat 
both  models  are  incorrect.  A  better  way  of  validating  this  range-dependent 

model  is  to  compare  its  solutions  to  experimental  measurements. 

11  12 

Ellis  and  Chapman,  ’  from  the  Defense  Research  Establishment  Atlantic 
(DREA),  Dartmouth,  Canada,  have  participated  in  a  sea-test  at  a  United  Kingdom 
continental  shelf.  One  of  the  test  areas  has  a  slight  range  dependence  of  the 
ocean  floor.  The  approximate  depth  of  the  bottom  is  100.0  mete’^s  and  the 
composition  of  the  bottom  is  mostly  chalk*  (see  Table  1)  witi.  a  few  meters  of 
sand  at  the  top.  They  modeled  this  ocean  environment  as  a  range-independent 
water  layer  over  a  semi-infinite  chalk  basement.  Inerefore,  they  neglected  the 
sand  sediment  and  the  ubiquitous  basalt  basement  that  should  be  located 
somewhere  under  the  chalk  sediment.  This  approximation  is  not  valid  at 
frequencies  much  lower  than  the  optimum  frequency  of  sound  propagation, 
because  the  penetration  capability  of  its  normal  modes  becomes  higher  and  they 
may  reach  the  depth  of  the  basalt  basement.  This  simple  two-layer  model  is 
also  not  valid  at  frequencies  higher  than  the  optimum  frequency  because  the 
effects  of  the  depth  dependent  water  column  becomes  of  paramount  importance  to 


*  Note  that  chalk  does  not  satisfy  the  second  inequality  of  Equations  (185), 
but  it  is  not  unusual  for  measured  attenuation  coefficients  to  be  highly 
erroneous  since  they  are  the  most  difficult  ones  to  obtain.  Even  though  these 
properties  of  chalk  may  by  questioned,  they  will  be  used  for  calculating  the 
transmission  loss  in  the  range-dependent  environment  described  in  References 
10  and  11. 
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the  transmission  loss  calculation.  They  found  the  optimum  frequency  for  the 
environment  to  be  in  the  vicinity  of  300.0  Hz.  At  much  higher  frequencies, 
the  effects  of  the  shear  waves  from  the  elastic  bottom  becomes  negligible 
compared  to  the  effects  of  the  depth-dependent  water  column.  Therefore,  we 
will  concentrate  on  the  frequencies  near  and  below  the  optimum  frequency. 

The  source  used  were  explosives  that  were  preset  to  detonate  at  a  depth 
of  37.5  ±  1  meters  and  a  hydrophone  was  located  at  71  meters  deep.  The  water 
depth  is  about  105  meters  at  the  location  of  the  hydrophone  and  it  has  a 
constant  slope  with  a  water  depth  of  95  meters  at  a  range  of  55  kilometers 
from  the  hydrophone.  This  corresponds  to  a  bottom  slope  of  0.01  degree  and  it 
can  be  considered  a  range- independent  wave  guide.  A  considerable  amount  of 
transmission  loss  measurements  have  been  provided  by  Chapman  for  the  1/3 
octave  band  center  frequencies  of  64,  128,  256,  512,  and  1024  Hz  as  a  function 
of  range  from  10  to  90  kilometers.  It  has  been  found  that  models  with  or 
without  shear  contribution  provide  nearly  the  same  transmission  loss  for 
frequencies  higher  than  256  Hz,  hence  the  lower  frequencies  will  be  considered 
here. 

Figure  31  is  the  three-dimensional  plot  of  the  range-dependent  coherent 
transmission  loss  for  a  frequency  of  128  Hz  and  a  source  depth  of  38  meters. 
The  bottom  is  a  semi-infinite  basement  of  chalk  with  the  properties  in  Table 
1.  The  water  column  has  a  sound  speed  of  1508  m/s  from  the  surface  to  a  depth 
of  28  meters,  and  a  constant  sound  speed  of  1494  m/s  from  a  depth  of  45  meters 
to  the  bottom.  The  density  is  a  constant  1  gm/cc  from  the  surface  to  the 
bottom.  The  range-dependent  wave  guide  has  been  divided  into  22  range- 
independent  range  segments.  The  first  few  depth  functions  are  T-R  modes 
similar  to  the  ones  for  a  sandy  bottom  (see  Figures  12  and  13)  and  all  but  the 
fundamental  mode  have  negligible  contribution  at  ranges  greater  than  about  20 
kilometers.  The  contour  plot  of  this  down-slope  wave  guide  is  shown  in  Figure 
32  where  the  slight  discontinuities  in  the  derivative  of  the  contours  are 
caused  by  the  range  segments. 

Comparison  of  the  measured  and  computed  transmission  loss,  for  the 

hydrophone  depth  of  71  meters,  is  displayed  in  Figure  33.  The  frequency  is 

128  Hz  and  the  range  is  extended  to  100  kilometers  to  accommodate  the  provided 

measurements.  This  transmission  loss  calculation  agrees  with  the  one  made  by 

11  12 

Ellis  and  Chapman  with  the  simple  two-layer  model.  ’  However,  their  model  is 
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overestimating  the  loss  at  frequencies  below  100  Hz. 

Figure  34  presents  the  calculated  and  measured  transmission  loss  for  the 
frequency  of  64  Hz.  Note  that  the  excited  normal  modes  are  rapidly  damped  at 
this  lower  frequency.  If  the  shear  speed  of  chalk  is  changed  to  the  fluid¬ 
like  value  in  Table  1,  the  calculated  coherent  and  incoherent  transmission 
loss  becomes  the  one  in  Figure  35.  Now  all  the  normal  modes  are  propagating 
with  much  less  attenuation,  but  this  fluid-like  model  is  underestimating  the 
loss.  Note  from  these  plots  the  tremendous  importance  of  the  shear  waves  in 
the  transmission  loss  calculation. 

However,  there  must  be  a  reason  for  the  disagreement  between  the 
theoretical  and  experimental  values.  Ellis  and  Chapman  speculated  that  a  deep 
reflector  is  causing  some  of  the  acoustic  energy  to  return  to  the  water 
column,  but  their  simple  model  is  incapable  of  including  more  layers. 

Under  the  assumption  that  their  suggested  deep  reflector  may  be  the 
omnipresent  basalt  basement,  a  semi-infinite  layer  has  been  included  in  our 
model  with  the  properties  of  basalt  given  in  Table  1.  The  depth  of  the  chalk- 
basalt  interface  was  taken  as  a  variable  in  order  to  fit  the  calculated 
transmission  loss  with  the  experimental  data.  This  inverse  scattering 
technique  provided  the  best  fit  for  a  chalk-basalt  interface  depth  of  about 
240  meters  and  the  resulting  transmission  loss  is  displayed  in  Figure  36.  The 
disagreement  at  ranges  greater  than  60  kilometers  may  be  due  to  the  extremely 
high  transmission  loss  that  makes  the  signal  fall  below  the  noise  level  of  the 
measured  data. 

Since  the  exact  location  of  the  sea-test  is  confidentially  kept  by  the 
Canadians,  there  is  no  way  we  can  verify  the  true  depth  of  the  basalt. 

However,  it  has  been  noticed  that  this  estimated  depth  is  typical  for  similar 
ocean  environments.  Also  it  has  been  verified  that  the  semi-infinite  layer  of 
basalt  hardly  changes  the  transmission  loss  calculations  at  128  Hz  (see  Figure 
33)  because  the  fundamental  depth  function  becomes  negligible  at  240  meters. 

It  is  true  that  the  wave  guide  just  considered  may  be  taken  as  range- 
independent.  Therefore,  an  up-slope  range-dependent  wave  guide  will  be  used  to 
test  our  model  for  a  steeper  bottom  slope.  A  25.0  Hz  source  is  located  at  a 
depth  of  112.0  meters  in  a  200.0  meters  deep  water  colum.n  of  constant  sound 
speed  (1500.0  m/s)  over  a  fluid-like  bottom  with  the  properties  in  Table  1, 
Three  trapped  and  four  radiating  modes  are  detected  in  this  range-independent 
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wave  guide.  The  source  has  been  placed  at  a  node  of  the  second  normal  mode  to 
avoid  its  excitation.  The  contour  plot  of  the  coherent  transmission  loss  is 
provided  in  Figure  37  and  its  three-dimensional  display  is  in  Figure  38  where 
all  seven  modes  have  been  included  in  their  calculation.  Since  the  Hankel 
function  computation  has  been  performed  using  the  asymptotic  approximation, 
the  near  field  (r  <  X  =  c/f  s  60  meters)  transmission  loss  is  not  correct  and 
is  not  displayed  in  the  given  plots. 

To  convert  this  wave  guide  into  a  range-dependent  one,  we  will  create  an 
up-slope  that  starts  at  five  kilometers  and  ends  at  ten  kilometers  from  the 
source  with  a  final  bottom  depth  of  150.0  meters  deep.  Beyond  ten  kilometers 
the  wave  guide  remains  range  independent.  The  slope  has  been  divided  into  50 
segments  and  the  third  trapped  mode  becomes  a  radiating  mode  in  the  shallow 
portion  of  the  wave  guide.  This  slope  has  an  angle  of  0.57  degree  and  only  the 
first  three  normal  modes  have  been  used  for  this  computation  because  the 
higher  order  modes  are  of  no  effect  to  the  transmission  loss  at  the  region  of 
interest.  The  results  in  this  range-dependent  wave  guide  are  given  in  Figures 
39  and  40.  Note  that  some  of  the  acoustic  energy  is  propagating  into  the 
bottom  as  a  consequence  of  the  slope  which  is  converting  the  third  trapped 
mode  into  a  radiating  mode  with  a  higher  imaginary  component  of  the 
eigenvalue.  Also  note,  by  comparison  of  the  range-dependent  case  (see  Figure 
39)  with  the  range- independent  case  (see  Figure  37),  that  the  transmission 
loss  near  the  source  is  almost  identical  to  the  one  for  the  range- independent 
case.  Hence,  assuring  the  proper  range-dependent  transmission  loss 
computations.  A  similar  propagation  behavior  was  detected  by  Jensen^ ^  and 
Miller^^ 

The  variation  of  the  real  part  of  the  three  eigenvalues  with  range 
segment  number  is  plotted  in  Figure  41  where  the  bottom  curve  with  the  highest 
variation  is  the  third  trapped  eigenvalue  as  it  becomes  a  radiating  one.  The 
imaginary  component  of  the  eigenvalue  is  displayed  in  Figure  42  where  the 
imaginary  part  of  the  third  mode  has  become  so  high  that  its  contribution  to 
the  transmission  loss  can  be  neglected.  The  real  part  of  the  third  normal  mode 
at  the  first  range  segment  with  the  water  depth  of  170  meters  is  given  in 
Figure  43.  As  the  water  depth  becomes  shallower  the  third  normal  mode  becomes 
the  radiating  mode  In  Figuie  44  foi  the  water  deptn  of  ISO  .neters.  At  150 
meters  water  depth  the  mode  develops  more  oscillations  into  the  bottom  (see 
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Figure  45)  and  its  amplitude  becomes  order  of  magnitude  smaller  displaying  its 
negligible  contribution.  Hence,  the  precipitated  propagation  of  the 
interference  pattern  in  Figure  39. 

So  far,  we  have  performed  transmission  loss  computations  for  range- 
dependent  wave  guides  with  negligible  shear  contribution  and  the  water  column 
has  been  a  single  layer  of  constant  sound  speed  and  density.  An  actual  sound 
speed  from  the  Arctic  Ocean  is  provided  in  Figure  2  with  its  salinity  and 
temperature  profiles.  The  very  low  temperature  of  the  environment  causes  the 
propagated  sound  to  be  much  slower  than  1500  m/s  and  the  high  fluctuations 
with  depth  are  caused  mainly  by  internal  currents  typical  of  the  shallow 
region  of  this  ocean.  This  microstructure  of  the  sound  speed  profile  is  highly 
important  to  the  propagation  of  high-frequency  sound.  For  a  25  Hz  source,  a 
valid  approximation  is  to  consider  the  200  meters  water  column  a  single  layer 
with  surface  sound  speed  of  1435  m/s  and  a  bottom  sound  speed  of  1460  m/s.  The 
sound  speed  gradient  is  0.125  s  ^  and  the  bottom  is  a  semi-infinite  layer  of 
sand.  To  consider  the  case  of  downslope  propagation,  the  bottom  depth 
increases  from  200  meters  at  five  kilometers  from  the  source  to  400  meters  at 
10  km  range.  The  resulting  transmission  loss  in  this  range-dependent 
environment,  with  the  bottom  slope  of  -2.29,  is  provided  in  Figures  46  and  47. 
Note,  by  comparison  with  Figures  16  and  17,  that  the  gradient  causes  the  sound 
to  interact  less  with  the  bottom,  therefore  causing  it  to  propagate  with  less 
loss.  Also,  as  the  bottom  depth  becomes  larger,  the  sound  gets  trapped  in  the 
surface  channel  caused  by  the  positive  gradient.  This  channeling  behavior  is 
also  modeled  by  ray  bundles  that  bend  upward  and  bounce  back  from  the 
pressure-release  surface  forming  caustics  at  the  regions  where  they 
contructively  interfere.  Finally,  note  the  destructive  interference  that 
occurs  in  the  bottom  at  about  seven  kilometers.  As  the  sound  bounces  from  the 
ocean  floor,  some  of  its  energy  gets  refracted  into  the  bottom.  However,  at 
the  range-dependent  region  of  the  wave  guide  the  angle  of  reflection  is 
affected  by  the  slope,  causing  most  of  the  reflected  energy  to  scatter  the 
surface  at  a  shallower  angle  and  become  trapped  in  the  water  column.  The 
combined  effects  of  sound  trapped  in  the  channel  and  the  reflections  from  the 
slope  contribute  to  the  easy  detection  of  surface  ships  and  submarines  from 
open-ocean  receivers. 

The  plots  for  the  range-dependent  wave  guides  do  exhibit  sound 
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penetration  into  the  bottom,  but  it  seems  to  be  more  dampen  than  the  one 

17  53 

reported  by  Jensen  and  the  one  by  Miller.  A  probable  explanation  is  that 
the  mode  coupling  terms  in  Equations  (240)  must  be  included  in  the  range- 
dependent  transmission  loss  calculations  since  the  adiabatic  approximation 
breaks  dovm  for  rapidly  varying  ocean  wave  guides. 

Another  important  step  for  a  better  ocean  model  is  to  include  the  effects 
of  axial  variations,  bince  the  number  of  radiating  modes  have  been  drastically 
reduced  with  the  approach  given  in  this  investigation,  the  problem  of  computer 
memory  and  storage  has  been  decreased  and  further  computations  can  be  pursued. 

This  work  represents  a  step  closer  to  the  final  three-dimensional  coupled 
normal-mode  model  with  shear  wave  from  the  ocean  floor  and  the  Arctic  snow/ice 
surface  layers. 
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CHAPTE3^  12 

CONCLUSIONS  AND  RECOMMENDATIONS 


A  new  sound  propagation  model  and  its  computer  code  has  been  developed 
based  on  the  theory  of  normal  modes.  This  normal  mode  model  has  been  expanded 
to  take  into  account  the  effects  of  the  elasticity  of  the  ocean  floor  and  the 
depth  dependence  of  the  acoustic  properties  by  dividing  the  wave  guide  into 
horizontal  layers  with  constant  density,  constant  shear  speed,  and  constant 
attenuation  coefficients.  However,  the  water  column  has  layers  of  linear  wave 
number  squared  to  better  simulate  the  sound  speed  profile.  It  has  been  found 
that  the  compressional  sound  speed  in  the  elastic  layers  can  also  have  linear 
wave  number  squared  and  the  density  in  the  water  layers  can  be  a  variable  and 
still  have  a  solvable  set  of  wave  equations.  However,  the  limited  knowledge  of 
the  detailed  properties  of  the  bottom  and  the  limited  applications  suggest 
that  these  flexibilities  can  be  excluded  from  the  computer  code.  Since  the 
absorptive  properties  of  the  bottom  is  so  high  and  the  attenuation  of  low- 
frequency  sound  in  the  water  is  so  low,  the  absorption  in  the  water  has  been 
neglected. 

The  newly  developed  normal  mode  model  searches  for  the  eigenvalues  in  the 
complex  wave  number  plane  using  the  Levenberg-Marquardt  algorithm  that 
searches  for  the  minima  of  the  magnitude  of  the  complex  determinant.  It  has 
been  found  that  the  absorptive  properties  of  the  semi-infinite  bottom  causes 
the  radiating  wave  number  spectrum  to  be  inherently  discrete,  hence  the  false 
boundary  introduced  by  Evans  has  been  eliminated  and  the  number  of  radiating 
modes  has  been  drastically  reduced.  The  reduced  number  of  modes  for  the 
transmission  loss  calculations  allows  for  the  feasibility  of  calculations  at 
higher  frequencies  and  deeper  ocean  wave  guides. 

The  transmission  loss  in  the  elastic  sediments  is  computed  using  the 
magnitude  of  the  acoustic  intensity  vector.  This  complex  intensity  vector  is 
the  scalar  product  of  the  stress  tensor  and  the  particle  velocity  vector.  The 
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intensity  vector  can  also  be  used  to  determine  the  direction  of  propagation  of 
the  acoustic  energy. 

Range  dependence  of  the  acoustic  properties  of  the  ocean  wave  guide  has 
been  taken  into  account  by  using  a  modified  version  of  the  adiabatic  normal 
mode  theory  to  include  the  shear  waves.  This  adiabatic  normal  mode  theory  has 
been  developed  with  the  assumption  of  a  slowly  varying  environment.  However, 
the  threshold  of  maximum  variation  is  not  known  because  of  the  high  degree  of 
complexity  of  this  multilayered  model  and  the  large  number  of  acoustic 
properties  that  can  be  varied  in  range. 

The  eigenvalues,  eigenfunctions,  and  range- independent  transmission  loss 
results  have  been  compared  to  the  benchmark  two-layer  model,  with  a 
semi-infinite  elastic  bottom,  by  Ellis  and  Chapman  and  the  perturbation  method 
for  fluid-like  bottom  by  Miller  yielding  excellent  agreement. 

The  range-dependent  coherent  transmission  loss  calculation  has  been 
compared  to  transmission  loss  measurements  by  the  Defence  Research 
Establishment  Atlantic  (DREA),  Dartmouth,  N.S.,  Canada,  at  the  United  Kingdom 
continental  shelf.  Very  good  agreement  was  obtained  at  128  Hz  and  above  with  a 
the  model  containing  a  semi-infinite  chalk  bottom.  However,  this  model 
overestimates  the  loss  at  frequencies  below  100  Hz.  At  these  lower  frequencies 
the  shear  and  compressional  depth  functions  extended  deeper  into  the  bottom 
where  acoustic  properties  are  unknown,  hence  the  ubiquitous  basalt  basement 
has  been  included  at  a  depth  of  240  meters  from  the  ocean  surface  to  provide 
the  agreement  at  64  Hz  without  changing  the  results  at  128  Hz  and  above. 

Hence,  this  multilayered  model  can  also  be  used  for  inverse  scattering 
purposes . 

The  up-slope  wedge-like  ocean  has  been  modeled  for  a  variable  slope  to 
observe  the  changes  in  the  transmission  properties  and  to  test  the  validity  of 
the  adiabatic  approximation.  Perfect  agreement  has  been  found  between  the 
range- independent  and  the  range-dependent  transmission  loss  when  all  the 
segments  had  the  same  acoustic  properties  and  layer  thickness.  As  the  slope 
increased  a  "tongue,"  similar  to  the  one  observed  by  Jensen  and  Miller,  was 
developed.  However,  increasing  inaccuracy  of  the  range-dependent  transmission 
loss  with  increasing  bottom  slope  is  expected  due  to  the  need  of  the  coupling 
terms  in  the  inhomogeneous  range  equation  which  involve  the  range  derivative 
of  the  eigenfunctions. 
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The  next  step  to  the  ultimate  transmission  loss  model  is  to  include  these 
coupling  contributions.  Also,  it  is  possible  to  include  the  azimuthal 
variations  of  the  ocean  using  the  adiabatic  approximation  to  obtain  a  three- 
dimensional  transmission  loss  model  for  slowly  varying  environments.  Finally, 
the  azimuthal  and  range  coupling  contributions  can  be  incorporated  in  the 
three-dimensional  model  to  simulate  sea-mounts  and  more  complex  underwater 
structures. 

Other  steps  to  improve  the  present  computer  code  are  to: 

1.  Include  the  effects  of  the  elastic  snow/ice  layers  at  the  surface  of  the 
ocean  model  to  simulate  the  sound  propagation  in  the  Arctic  environment. 
With  such  a  model,  it  is  possible  to  study  the  effects  of  shear  waves  on 
ice-mounted  receivers. 

2.  Include  layers  with  linear  variation  of  the  density  with  depth.  The 
variation  of  density  with  depth  in  the  water  column  has  been  measured  and 
found  to  be  of  minimum  importance,  but  its  variation  in  the  elastic 
sediments  is  often  of  considerable  importance. 

3.  Include  elastic  layers  with  variable  compressional  sound  speed.  It  has 
been  theoretically  proven  in  this  work  that  the  elastic  wavs  equation 
representing  a  layer  with  variable  compressional  sound  speed  can  be 
separated  into  an  equation  for  shear  waves  and  one  for  the  compressional 
waves . 

4.  Include  absorption  effects  from  the  water  column.  The  absorption  in  the 
water  at  low  frequencies  is  negligible.  However,  its  contribution  at 
higher  frequencies  becomes  important  and  it  must  be  included  in  the 
transmission  loss  calculation. 

5.  Include  the  effects  from  surface  and  bottom  roughness.  Very  simple 

7  4 

equations  have  been  derived  by  Kuperman  and  Ingenito  with  the  Kirchhoff 
approximation.  The  equations  ignore  the  contributions  from  the 
non-specularly  reflected  acoustic  energy  and  they  may  be  added  to  the 
imaginary  part  of  the  complex  wave  numbers  after  the  eigenfunctions  are 
computed. 

6.  Refine  the  searching  algorithm  to  "guarantee"  the  uniform  convergence  to 
all  the  complex  eigenvalues.  The  present  searching  algorithm  may  not  be 
able  to  find  all  the  complex  eigenvalues  for  a  water  column  with  two  or 
more  channels  since  these  create  degenerate  eigenvalues  with  irregular 
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spacings. 

7.  Include  an  option  to  obtain  the  transmission  loss  as  a  function  of 

frequency  and  to  account  for  the  frequency  spectrum  of  the  signal  emitted 
by  the  source  (its  signature).  The  current  computer  code  calculates  the 
transmission  loss  of  continuous  wave  (CW)  acoustic  signals  and  extra 
computations  are  required  to  obtain  the  transmission  loss  of  pulses  and 
other  wave  forms. 
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FIGURE  1 


F(x)=Ap(x) 


<-  dx  ^ 


F ( x+dx ) =Ap ( x+dx ) 


area 


b. 


V  (x)  dt 

X 


dx  -> 


A  = 


V  (x+dx)  dt 
area  x 


MODELS  TO  D’T^IVE  THE  EULER  EQUATION  OF  MOTION  (a)  AND  THE 
CONTINUITY  EQUATION  (b) 
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1420.  □  SOUND  SPEED  CM/S3  1480. 
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FIGURE  2.  SOUND  SPEED,  SALINITY,  AND  TEMPERATURE  PRCFILES  TAKES  IN  THE  EAST 
GREO^LAf^D  C'.T’RENT  (78°N) 
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FIGURE  10.  THE  COHERENT  (SOLID  CURVE)  AND  INCOHERENT  (DASH  CURVE) 

transmission  loss  versus  range  for  a  25  HZ  SOURCE  IN  A  WATER 
COLUMN  200  METERS  DEEP  OVER  A  FLUID-LIKE  SEMI-INFINITE  BASEM 
(THE  SOURCE  AND  RECEIVER  DEPTH  IS  112  METERS) 


PEAL  PART  nf  r-lODE  tt 
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figure  16.  A  THREE-DIMENSIONAL  PLOT  OF  THE  COHERENT  TRANSMISSION  LOSS  (DB)  AS 
A  FUNCTION  OF  RANGE  (KM)  AND  DEPTH  (M)  FOR  A  25  HZ  SOURCE  AT  112  M 
DEPTH  IN  A  200  M  WATER  COLUMN  OVER  A  SEMI -INFINITE  SAND  BASEMENT 
(ALL  SEVEN  MODES  WERE  USED) 
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FIGURE  23.  LOCATION  OF  THE  COHPLEX  EIGENVALUES  FOR  A  25  HZ  SOURCE  IN  A  WATER 
COLUMN  200  METERS  DEEP  OVER  A  SEMI-INFINITE  BASALT  BASEMENT 
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FIGURE  28.  THE  REAL  PART  OF  THE  LAST  MODE  FOR  A  25  HZ  SOURCE  IN  A  WATER 
COLUMN  200  METERS  DEEP  OVER  A  SEMI-INFINITE  BASALT  BASEMENT  C 
SOLID  CURVE  IS  THE  OOMPRESSIONAL  EIGENFUNCTION  AND  THE  DASHED 
CURVE  IS  THE  SHEAR  EIGENFUNCTION) 
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A  FUNCTION  OF  RANGE  (KM)  AND  DEPTH  (M)  FOR  A  128  HZ  SOURCE  AT  38  M 
DEPTH  IN  A  SLIGHTLY  RANGE-DEPENDENT  (DOWN-SLOPE)  WATER  COLUMN  OVER 
A  SEMI- INFINITE  CHALK  BASEMENT  (TEN  MODES  WERE  USED) 


FIGURE  33.  COMPARISON  WITH  EXPERIMENTAL  MEASUREMENTS  FOR  A  128  HZ  SOURCE  IN  A 
SLIGHTLY  RANGE-DEPENDENT  (DOWN-SLOPE)  WATER  COLUMN  OVER  A 
SEMI-INFINITE  CHALK  BASEMENT 


FIGURE  34.  COMPARISON  WITH  EXPERIMENTAL  MEASUREMENTS  FOR  A  64  HZ  SOURCE  IN  A 
SLIGHTLY  RANGE-DEPENDENT  (DOWN-SLOPE)  WATER  COLUMN  OVER  A 
SEMI -INFINITE  CHALK  BASEMENT 


FIGURE  35.  COMPARISON  WITH  EXPERIMENTAL  MEASUREMENTS  FOR  A  64  HZ  SOURCE  IN  A 
SLIGHTLY  RANGE-DEPENDENT  (DOWN-SLOPE)  WATER  COLUMN  OVER  A  SEMI- 
INFINITE  "CHALK"  BASEMENT  WITH  A  NEGLIGIBLE  SHEAR  SPEED  (THE  SOLID 
CURVE  IS  THE  COHERENT  TRANSMISSION  LOSS  AND  THE  DASHED  CURVE  IS 
THE  INCOHERENT  TRANSMISSION  LOSS) 
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FIGURE  39.  A  CONTOUR  PLOT  OF  THE  COHERENT  TRANSMISSION  LOSS  (DB)  AS  A 

FUNCTION  OF  RANGE  (KM)  AND  DEPTH  (M)  FOR  A  25  HZ  SOURCE  AT  112  M 
DEPTH  IN  AN  UP-SLOPE  WATER  C01.UMN  WITH  A  0.57  DEGREE  SLOPE  OVER  A 
SEMI -INFINITE  FLUID-LIKE  BASEMENT  (THE  FIRST  THREE  MODES  WERE  USED 
AND  THE  WAVE  GUIDE  WAS  DIVIDED  INTO  50  SEGMENTS) 
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FIGURE  41.  THE  REAL  COMPONENT  OF  THE  FIRST  THREE  EIGENVALUES  AS  A  FUNCTION  OF 
THE  RANGE  SEGMENT  NUMBER  FOR  A  25  HZ  SOURCE  IN  AN  UP-SLOPE  WATER 
COLUMN  WITH  A  0.57  DEGREE  SLOPE  OVER  A  SEMI-INFINITE  FLUID-LIKE 


FIGURE  43,  THE  REAL  PART  OF  THE  THIRD  DEPTH  FUNCTION  FOR  A  25  HZ  SOURCE  IN  A 
170  M  WATER  COLUMN  OVER  A  FLUID-LIKE  SEMI-INFINITE  BASEMENT 
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FIGURE  46.  A  CONTOUR  PLOT  OF  THE  COHERENT  TRANSMISSION  LOSS  (DB)  AS  A 

FUNCTION  OF  RANGE  (KM)  AND  DEPTH  (M)  FOR  A  25  HZ  SOURCE  AT  112  M 
DEPTH  IN  A  WATER  COLUMN  WITH  SOUND  SPEED  GRADIENT  OF  0. 125  S'^  AND 
A  BOTTOM  SLOPE  OF  -2.29  OVER  A  SEMI-INFINITE  SAND  BASEMENT  (SEVEN 
MODES  WERE  USED  AND  THE  WAVE  GUIDE  WAS  DIVIDED  INTO  50  SEGMENTS) 
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FIGURE  47.  A  THREE-DIMENSIONAL  PLOT  OF  THE  COHERENT  TRANSMISSION  LOSS  (DB)  AS 
A  FUNCTION  OF  RANGE  (KM)  AND  DEPTH  (M)  FOR  A  25  HZ  SOURCE  AT  112  M 
DEPTH  IN  A  WATER  COLUMN  WITH  SOUND  SPEED  GRADIENT  OF  0.125  S'*  AND 
A  BOTTOM  SLOPE  OF  -2.29  OVER  A  SEMI-INFINITE  SAND  BASEMENT  (SEVEN 
MODES  WERE  USED  AND  THE  WAVE  GUIDE  WAS  DIVIDED  INTO  50  SEGMENTS) 
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TABLE  1.  GEO- ACOUSTIC  PROPERTIES  OF  THE  VARIOUS  SEDIMENTS 


Bottom 

Type 

Density 

Compressional 

speed 

Compressional 

attenuation 

Shear 

speed 

Shear 

attenuation 

(gm/cc) 

(m/s) 

(dB/kHz-m) 

(m/s) 

(dB/kHz-m) 

Fluid-like 

1.15 

1704.5 

0.29 

1.0 

1.00 

Clay-Silt 

1.60 

1515.0 

0.50 

100.  0 

1.00 

Sand 

2.00 

1800.0 

0.70 

600.0 

1.50 

Basalt 

2.60 

5250.0 

0.20 

2500.0 

0.50 

Chalk 

2.20 

3200.0 

0.  10 

1000.0 

1.00 
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TABLE  2 


COMPARISON  OF  THE  CALCULATED  REAL  PART  OF  THE  TRAPPED  EIGENVALUES 


Mode# 

Rigid  Model 

Soft  Model 

Perturbation 

Exact  Model 

1 

0. 1044248 

0. 1035350 

0. 1041654 

0. 1040583 

2 

0. 1020346 

0.0998963 

0. 1012651 

0. 1014331 

3 

0.0970778 

0.0935177 

0.0963706 

0.0960692 

4 

0.0891272 

0 . 0837758 

5 

0.0772641 

0.0692656 

6 

0.0591806 

0.0456463 

7 

_ 

0.0232692 
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